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Abstract 

The objective of this thesis is to further develop and 
test a stochastic model of turbulent combustion in 
recirculating flows. There is a requirement to increase 
the accuracy of multi-dimensional combustion predictions. 

As turbulence affects reaction rates, this interaction must 
be more accurately evaluated. In this work a more 
physically correct way of handling the interaction of 
turbulence on combustion is further developed and tested. 

As turbulence involves randomness, stochastic modeling is 
used. Averaged values such as temperature and species 
concentration are found by integrating the probability 
density function (pdf) over the range of the scalar. The 
model in this work does not assume the pdf type, but solves 
for the evolution of the pdf using the Monte Carlo solution 
technique. The model is further developed by including a 
more robust reaction solver, by using accurate 
thermodynamics and by more accurate transport of elements. 
The stochastic method is used with Semi-Implicit Method for 
Pressure-Linked Equations. The SIMPLE method is used to 
solve for velocity, pressure, turbulent kinetic energy and 
dissipation. The pdf solver solves for temperature and 
species concentration. Thus, the method is partially 
familiar to combustor engineers. The method is compared to 
benchmark experimental data and baseline calculations . The 
baseline method was tested on isothermal flows, evaporating 
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sprays and combusting sprays. Pdf and baseline predictions 
were performed for three diffusion flames and one premixed 
flame. The pdf method predicted lower combustion rates than 
the baseline method in agreement with the data, except for 
the premixed flame. The baseline and stochastic predictions 
bounded the experimental data for the premixed flame. The 
use of a continuous mixing model or relax to mean mixing 
model had little effect on the prediction of average 
temperature. Two grids were used in a hydrogen diffusion 
flame simulation. Grid density didn't affect the 
predictions except for peak temperature and tangential 
velocity. The hybrid pdf method did take longer and 
required more memory, but has a theoretical basis to extend 
to many reaction steps which cannot be said of current 
turbulent combustion models. 
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Chapter 1 
INTRODUCTION 
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1.1 Motivation for the thesis 

Improving the numerical simulation of combustor flows 
is a never-ending endeavor. There is a continuing demand to 
increase gas turbine performance and decrease emissions. 

These demands do not complement each other. The demand for 
increased performance is typically met by increasing cycle 
pressure and final temperature. This is allowed by 
continual advances in metallurgy. However , increasing 
pressure and temperature usually increases emissions. Thus, 
combustor designs are being evaluated to test for increased 
performance with lower emissions. One method to help in 
designing gas turbine combustors is Computational Fluid 
Dynamics and combustor models. These models are benchmarked 
or validated using existing engine or test data. CFD 
combustor model predictions are accurate or calibrated _ 
within similar engine designs. Thus, new concepts cannot 
be totally evaluated by current CFD models. The number of 
prototypes in an engine development program may be reduced 
by using CFD combustor model predictions. Initial 
computational designs must be tested on engine stands. CFD 
can be also be used to correct design deficiencies found in 
prototypes. Commercial gas turbines are very difficult and 
costly to instrument and test. The use of CFD models 
provides data at all locations in a combustor field. Thus, 
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the nature of CFD complements engine test programs. 

Newer engine designs are going to higher peak 
temperature and pressure to improve both fuel economy and * 

thrust to weight ratio. Severe reductions in emissions are 
being required by governments . Some European countries are 
considering instituting fees based upon the total amount of 
pollutants emitted by aircraft flying through their 
airspace. Older aircraft producing more pollution would be 
more heavily taxed. The evaluation of proposed engine 
designs is causing a demand for increased accuracy in CFD 
combustor models as well as a wider range of validity in 
those models. Engine designers are tolerating more complex 
combustor models in a desire for more accurate CFD 
predictions. Computer resources to allow more accurate 
modeling is becoming more affordable with time. Combustion 
simulations were once entirely done on mainframe computers 
or supercomputers. Running many combustion simulations on 
supercomputers can cost thousands of dollars. Workstations 
have been developed to the point that taking into account 
the amount of time jobs are in queues, the total turn-around 
time is comparable for high-end workstations and * 

supercomputers. The cost of doing calculations on 

f 

workstation is significantly less than running on 
supercomputers. To further reduce turn-around time of 
workstations, computers may be clustered together to give 
performance similar to dedicated super-computers. Many 
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workstations see relatively little use at off hours. With 
suitable networking, these under-utilized workstations may 
be used to perform complex 3-D predictions at small 
additional cost. 

— Progress in the Simulation of Combustion 
Much of the model improvement in the last decade has 
been in the area of gridding. Curved boundaries and round 
dilution jets are modeled without stair-stepping the 
geometry. Unfortunately, generating grids for complex 
combustors can take up most of the time to perform the 
calculation. One way around this is to overlap simpler 
grids. This technique essentially allows imbedding grid for 
a flow feature into another simpler grid which poorly 
defines the flow feature. A good example is a cylindrical 
grid for a rod imbedded into a cartesian grid for a 
rectangular duct. Another technique being developed is the 
use of unstructured grid solvers. In unstructured grids, 
the grid is generated by breaking the flow geometry into 
many smaller general elements. If more refinement is 
needed, the elements are further divided. The data storage 
for fluid elements is quite different than for structured 
grid solvers. The unstructured grid method incurs 
additional work to keep track of neighboring grid cells. 
This technique is only in its infancy. It is unknown which 
solution method will prove best for combustor design. 

Over the past twenty years, most elliptic flow 



calculations have been done using variants of the Semi- 
Implicit-Method for Pressure-Linked-Equations developed by 
Patankar 1 . The usual implementation is to assume low Mach 
number flow. In this flow regime, pressure is largely a 
relative variable, that is, pressure is fairly constant. 
Changes in density are due to changes in molecular weight 
and temperature. This solution technique is referred to as 
pressure based method. Industry regards the solution method 
as robust, but, as rather dissipative. If flow is at an 
angle to the grid, calculated flow quantities may be 
dissipated. A great deal of effort was expended in the 
eighties to develop improved differencing methods to improve 
accuracy of pressure based schemes. Skew differencing was 
developed to take into account the angle of the flow to the 
grid. Quadratic upwind interpolation and second-order 
upwind differencing were also developed. Theoretically 
these schemes are more accurate, but, the schemes are not as 
robust as first-order upwind differencing schemes. 
Convergence is much more difficult for improved schemes. 

Some of these schemes blend in varying amounts of upwind 
differencing to improve robustness. 

Density based solvers are thought by some individuals 
to be less dissipative than the pressure based solvers. 
However, these solvers must include artificial dissipation 
terms for robustness. Density based solvers are efficient 
in the solution of higher Mach number or compressible flows. 


5 


The convergence of the solver greatly deteriorates with 
decreasing Mach number as the flow becomes incompressible. 
Solution techniques using density based solvers for low 
speed flows are being developed, such as NASA Lewis Research 
Center's ALLSPD 2 code. The ALLSPD code uses 

preconditioning, gauge pressure and pseudo- time stepping to 
maintain good convergence rates at low Mach number flows. 

This solution technique is in the development phase. 

Another technique to improve convergence rates is to 
use multi-grid solvers. Joshi and Vanka 3 performed multi- 
grid calculations for gas-turbine type geometries. The 
method was also demonstrated for 3-D geometries for studying 
hot gas ingestion for short take off and vertical landing 
aircraft 4 . Multi-grid was also used in a numerical ramjet 
study by Vanka and Krazinski 5 . A four-step combustion 
mechanism with combination eddy-breakup and arrhenius 
reaction rate submodel was used in the study. A single grid 
system was used in the solution of species concentration due 
to the strong influence of the source terms in the chemical 
kinetic equations. 

The overall chemistry or reaction kinetics of fuels can 
be very complicated. Also, combustion in practical gas 
turbine combustors is highly turbulent. Improved turbulence 
models for variable density flows are the least developed of 
all turbulence models, and typically have been calibrated 
for isothermal incompressible flows. Turbulence is usually 



ignored by researchers who study fundamental combustion 
processes the most, the experts in chemical reaction 
kinetics. These experts predict reaction rates accurately 
over the range of combustor temperature, pressure and 
species concentration. Unfortunately, this is accomplished 
by employing dozens of chemical reactions with widely 
varying time characteristics. These kinetic reaction 
systems must be solved by specialized solvers using small 
time steps, predictor-corrector steps, etc. Many reactions 
involve species the engineer isn't concerned with. Not 
enough effort is being spent on producing reduced chemical 
reaction mechanisms that will give acceptable engineering 
predictions. Laminar combustion experts have yet to develop 
simpler kinetic reaction schemes which they feel are 
accurate enough to predict exhaust emissions for different 
engine concepts. Typical combustor models employ five to 
ten species and as few chemical reactions as possible. The 
chemical kineticists typically insist on twenty or more 
reactions for hydrocarbons and would like to see as many as 
fifty reactions. The time steps needed to accurately solve 
species concentrations for these kinetic mechanisms can be 
extremely small, much shorter than the pseudo time step 
taken in typical combustor calculations. A good time step 
for the LSENS 6 general kinetics and sensitivity solver is 
of the order of ten micro-seconds. The code employs 
variable time steps to improve solution speed. The solution 
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efficiency has been substantially improved from the original 
solver. This code is one of the most highly developed and 
efficient chemical kinetics solver existing today. However, 
it is still too computationally expensive to use for solving 
complex multi -dimensional combustor flows. Techniques used 
to accurately solve laminar one dimensional combustion are 
too complex and costly at the present time to use for the 
solution of multi-dimensional turbulent recirculating flows. 
Typical design models use a very small number of reaction 
steps to predict performance. Engine emissions are usually 
predicted by post-processing, adjusting constants or 
employing empirical corrections. 

1.3 - Turbulent Combustion 

The prediction of combustion is greatly affected by 
turbulence. The length of current gas turbine combustors is 
fairly short. Fuel and oxidant must be thoroughly mixed 
within a short distance at fairly high velocities. 

Combustor velocities are greater than laminar flame 
velocities. Combustor exit velocities in the newer gas 
turbine engines are getting into the high subsonic realm. 
The overall mixing of fuel and oxidants in a combustor is 

greatly augmented by turbulence. 

Turbulence causes changes in the instantaneous values 
of velocity, species and temperature. Reaction rates are 
nonlinear functions of concentration and temperature. Thus 
the averaged reaction rate isn't the same as the reaction 



rate of averaged species concentration and temperature. 
Previous modeling has employed many ad-hoc corrections to 
account for the interaction of turbulence on combustion. 

One such correction is to include a factor to convert 
laminar reaction rates into turbulent reaction rates. 
Turbulence produces eddies or modules with varying 
concentrations of fuel and oxidant. These eddies go through 
a breakup or decay process before the eddies completely mix. 
As turbulence continually generates new eddies in the flows, 
this is an ongoing process. Thus, turbulence causes 
incomplete mixing which causes an overlap of fuel and 
oxidant species concentration particularly near 
stoichiometric conditions where there is the greatest 
chemical reactivity. The overlap in average reactant specie 
concentration makes turbulent combustion resemble finite 
rate combustion. Current turbulent combustion methods do a 
good job of predicting combustor exit temperature patterns 
and major species concentrations. Minor species and 
intermediates, especially carbon monoxide, are not well 
predicted. Unfortunately, minor species can cause 
significant pollution levels. The level of engine emissions 
is being monitored and future emissions have been regulated 
due to serious environmental considerations. The 
development of an American supersonic transport was canceled 
in the seventies due to concerns of immense nitrous oxide 
production at high altitudes. Lakes and forests are being 
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adversely affected by current fossil fuel combustion. The 
combustion process must be better predicted to aid in the 
design of lower emissions combustors. 

1.4 - Stochastic Turbulent Combustion Models 

A more physically correct model of the turbulence- 
chemistry interaction is to treat combustion as a stochastic 
process. This is the idea behind current assumed shape 
probability density function models. These models assume 
certain shapes for the probability of a conserved scalar. 

The probability density function may be inferred by solving 
modeled equations for the means and variances of variables . 
The major assumption in presumed shape pdf modeling is that 
reaction rates are fast in these models. This can be 
partially relaxed by including another independent variable, 
typically for reactedness. The additional variable is 
usually assumed to be uncorrelated to the conserved scalar. 
Unfortunately, the additional variable is usually 
correlated. Including more species and chemical reactions 
adds complexity and increases the number of assumptions. 

Even allowing for incomplete or partial combustion in 
assumed shape pdf models causes predictions of high reaction 
completeness as the simplified combustion mechanisms 
employed typically used are only valid at high temperatures. 

Another stochastic method involves solving for the 
evolution of probability density function. This modeling 
eliminates the chemical reaction rate closure problem 



inherent in traditional combustion modeling. Stochastic 
treatments of species transport are said to exactly include 
the chemical reaction rate term. This has been known for 
quite some time, but, the capability to perform stochastic 
calculations of practical combustor flowfields is being 
developed at the experimental or research level. Actual 
engine combustor design predictions have to be done in a 
short period of time to be effective in the combustor design 
process. Also, industrial predictions have to be cost 
effective. Stochastic modeling demands large expensive 
computer resources . As computer power continues to increase 
and the cost of doing combustion simulations becomes more 
acceptable, the use of stochastic combustion models will 
increase . 

1.5 - Contribution of the Present Work 

Assumed shape pdf models have been used for over a 
decade and development of these models has largely 
stagnated. The next step in pdf modeling is the composition 
pdf, which requires roughly an order of magnitude increase 
in memory and computer time. Pdf modeling allows a 
substantially more robust and accurate method of calculating 
minor species for turbulent combustion, albeit at 
substantially increased cost. 

Newly developed hybrid pdf models have predicted 
excellent agreement with experimental data for parabolic or 
one-way combustor flow. The method is called hybrid as 
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stochastic methods are used to solve for average species 
composition, density and temperature, while traditional CFD 
techniques are used to solve for other flow quantities. 

Full stochastic simulations of all flow quantities is beyond 
current capabilities. The hybrid method retains the 
traditional CFD techniques which engineers are used to, 
while it treats turbulent chemistry interactions in a much 
more physically correct manner . 

The purpose of my work is to incorporate and further 
develop a composition and enthalpy pdf model into a 
computational fluid dynamics model. The hybrid pdf model 
uses individual species particles to solve for the evolution 
of composition and temperature pdfs. The hybrid pdf model 
consists of consecutively solving for velocity, pressure, 
turbulence length and time scales, species transport and 
chemical reaction. A Monte Carlo solution technique is used 
to predict species transport, mixing and species reactions. 
Many improvements to the original pdf model were made in the 
work. The pdf model was improved to more properly calculate 
species transport. The original model used an incorrectly 
coded upwind differencing scheme to calculate species 
transport. A more accurate scheme was implemented. The 
original model had major thermodynamic simplifications. 

This was corrected to better predict actual combustor data. 
The original model had an extremely simple combustion 
solver. This solver worked with the original reaction rate 
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constants, but, didn't for the reaction models in this work. 
An improved species solver was implemented. Pdf results are 
properly compared against good CFD results. Some pdf 
predictions are performed against CFD predictions done years 
ago, and claiming the improved results are due solely to pdf 
modeling. This ignores many improvements which have been 
made to CFD models over the years. Pdf predictions have 
also been compared to initial or other poorly refined 
predictions. The hybrid pdf model will be tested against a 
proven baseline turbulent combustion model and experimental 
data for various research combustor flowfields. The 
baseline model uses the pressure based method to solve for 
velocity, pressure, mean species, species fluctuation and 
turbulence scales in elliptic flow fields. The k-e 
turbulence model is used. A fuller description of the model 
development in this work is covered in chapter 4. 


? 
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Chapter 2 

BASELINE SOLUTION METHOD FOR CONVECTION-DIFFUSION PROBLEMS 
The Semi-Implicit Method for Pressure Linked Equations 
or SIMPLE method of Patankar is used in conjunction with the 
hybrid composition pdf method in this thesis. The SIMPLE 
algorithm method is used for solving momentum, overall mass 
conservation, turbulent kinetic energy and dissipation rate 
of turbulent kinetic energy. The pdf method used here 
solves turbulent combustion and species transport. This 
chapter covers a rather short discussion on the SIMPLE 
method, boundary conditions used and source terms for 
combustion species. 

2.1 - Governing Equations 

The gas phase equations can be put in the following 
general form for cartesian geometries: 




■p ffi) 


+ 



T* d<t> 


= 5 * 


( 2 . 1 ) 


where p is the density, U is the axial velocity, V is the 
radial velocity, is the appropriate diffusion coefficient 
and S^, is the source term for the variable <(>. For example, 
the diffusion coefficient for the laminar momentum equations 
would be p 2 , the laminar viscosity. Turbulence introduces 
additional complexity into the transport equation. The way 
this is traditionally done for incompressible flow is to 
perform a Reynolds decomposition of velocity into a mean and 
fluctuating components: 
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U = V + u' 


where, 


( 2 . 2 ) 


£7 = 



t 0 +dt 

J Udt 

to 


( 2 . 3 ) 


Substituting in the Reynolds decomposition of axial and 
radial velocity into the momentum equation and time 
averaging results in terms such as: 

pTFiP, pvV, and puV', etc 

which are known as Reynolds stresses. The first two terms 
are normal stresses. These can be incorporated into the 
pressure term. If density varies, the dilatation or 
divergence of velocity isn't zero. The dilatation can also 
be included in the pressure term: 

p * i 1 ^) * 12 - 4) 

where |i c is the turbulent coefficient of viscosity. The 
Reynolds stresses are unknown. The Reynolds stresses must 
involve modeling at some point due to problems with closure. 
Lower level modeling involves simple correlations for the 
stresses. Higher level or second moment modeling involves 
solving modeled equations for the transport of Reynolds 
stresses. These equations are formed by a sequence of 
multiplication of the momentum equations, Reynolds 
decomposition and time averaging. These equations have 
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unknown triple correlation terms , which must then be 
modeled. Lower level turbulence modeling has proved 
suitable for prediction of turbulent combustion. The non- 
normal Reynolds stresses are modeled using a generalization 
of Boussinesq's eddy viscosity concept or the gradient model 
assumption : 


■ij 




dul 

dx* 


du^ 


dx 


i/ 




( 2 . 5 ) 


Following the Prandtl -Kolmogorov definition of eddy 
viscosity, the viscosity coefficient is composed of length 
and velocity scales. The preferred method of supplying 
these scales is the two equation k-e turbulence model. The 
k-t model was developed by Jones and Launder 7 . The term k 
is known as the turbulent kinetic energy: 


k = 



( 2 . 6 ) 


The term 8 is known as the turbulent kinetic dissipation 
rate. If the characteristic length is defined as k 3/2 /e and 
velocity as k 1/2 , the turbulent diffusion coefficient can be 
defined as: 


H C = 



( 2 . 7 ) 


Equations can be formed for k and £. However, these 
equations have been modified to better predict turbulence. 
The modeled equations are: 
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tc ,p »*4 


puk-r k 


dk 

3x 


dy 


pVk-T k 


dk 

dy 


= - pe 


( 2 . 8 ) 


and 

~§t (p c > * H pC,e - r --i 

-|^[pve-r e |l] = C^Gt/k - C 2 pe 2 /k 
where G is the turbulence generation term: 


( 2 . 9 ) 


dUi 

G ’ 


( 2 . 10 ) 


The k-e model allows for the variability of length and 
velocity scales with history effects. That is, the 
viscosity isn't just a function of local gradients or 
scales . 

Some of the non-normal turbulent stresses can be 
incorporated into diffusion terms in the general equation by 
defining a effective diffusion coefficient which is the sum 
of the laminar and turbulent diffusion coefficients. 

= Pj + Pc (2.11) 

The left over turbulent terms are put into source terms for 
the momentum equations. For example, the axial momentum 
source term for Cartesian flows is: 



( 2 . 12 ) 
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2.2 - Discretization of the general equation 

Except for simplified problems, the general equations 
do not have an analytic solutions. The equations must be 


solved numerically. The 
equation is discretized by 
using the control volume 
approach. The general 
transport equation is 
integrated over the dashed 
control volume for the node 
P shown in figure 2.1. The 
values at the faces of the 
control volume are presumed 
to preside over the whole 
face. This gives: 



figure 2.1 Control Volume 


Fjb w +D W {Q W -Q p ) +f At>- d n^ n -$p) ~ 

f b * b +d s (* p -* s ) - 


(2-13) 


= (Sc+S^p) axa y 


where the source term has been linearized and the convection 
coefficients at the control volume faces are given by: 


F e = (pm o Ay, 


(2.14a) 


F n = (p 


(2.14b) 


F a = (pV) D AX, 


(2.14c) 
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F g = (pV) gAX, (2 . 14d) 

where Ax and Ay are the width and height of the control 
volume. The diffusion coefficients are given by: 





(2.15a) 


= 



(2.15b) 


(2.15c) 







(2 . 15d) 


where the term (Ax) e is the distance between the points P 
and E, etc. 

The equation is further developed by manipulation of 
the continuity equation. Integrating the continuity 
equation over the same volume gives: 

(pp-pf) + p e _ F w + f d _ p s = o (2.16) 

A t 

If the integrated continuity equation is multiplied by <J) P 
and subtracted from the integrated general equation, the 
following equation is obtained: 
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/a j.o \ pP°AXAy (2.17) 

(<Pp -< Pp) — — 

+F g (<|> e -<j>p) -D E (b E -4> P ) -F w ($ w -b P ) +d k W p - 4> w ) 

+F a ($ D -$ P ) -D n ($ n -$ p ) ~F b (<|> s -<l>p) +D S {$ P -$ S ) 

= {S C +S^p) AXAy 

This equation involves intermediate values at the cell 
faces. The goal is to discretize the above equation into an 
equation involving <J> P and its neighboring nodes: 

a p <t> p = a& B + a$„ + a^ N + a£ s + b (2.18) 

This may be done by forming interpolations involving the 
adjoining nodes for the intermediate values at the cell 
faces. However, if any of the above coefficients is allowed 
to become negative, the answers may become non-physical. To 
ensure physically realistic results the hybrid differencing 
scheme is used in this work. The hybrid differencing scheme 
is a combination of central differencing and upwind 
differencing. Hybrid differencing uses more accurate 
central differencing when the absolute value of the ratio of 
the convective flux to the diffusive flux (or Peclet number) 
is less than two. When the Peclet number is greater than 
two, the upwind differencing scheme coefficients are used, 
which is physically realistic at high convection rates. 

The use of hybrid differencing is illustrated with the 
following expression for fluxes on the east side of the 


control volume: 



20 


F e <4> e -<!>?) -d e (4>£r-4> P ) = a £ (<|> p -^ £ ) (2.19) 

If the absolute value of the Peclet number is less than two, 
central differencing is used for the term <)> e , thus: 



D e ($> E <J)p) — a £ (<|)p <|)p) 


( 2 . 20 ) 


which gives : 


a E = D e --^, for |Pe|*2 (2.21) 

If the Peclet number is greater than two, the convective 
flux from node P overwhelms the diffusive flux. The 
diffusive flux is ignored and <J> e is taken as <t> P , and: 


F e (<t>p-<i>p) = a, ($*-**) 


( 2 . 22 ) 


which gives : 

a E = 0, for Pe > 2 (2.23) 

If Pe is less than -2, then the convective flux from node E 
is much larger than the diffusive flux, which is ignored. 
Thus : 


which gives : 


^ e (4> £ -<l>p) = a £ (<(>p-4> B ) 


(2.24) 


a E = -F e , for Pe < -2. (2.25) 

The general value for a E has been worked out to be: 


= 





0 


(2.26) 
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where the symbols I ] stands for the largest of the 
quantities within it. The formula for a w is: 


A w 



(2.27) 


The formulas for a N and a s are similar. The formula for the 
rest of the coefficients are: 


a p = a E + a w + a N + a s + ap - Sp&x&y, 


(2.28a) 


ppAXAy 
At ' 


(2.28b) 


b = ScAXAy + ap4>p (2.28c) 

2.3 - Solution for Pressure 

An offset or staggered grid is used to solve the 
momentum equations. The nodes for axial velocity occur on 
the east and west faces of the scalar control volume. Thus, 
the control volume for axial momentum is offset from the 
■scalar control volume. Thus the discretization equation for 
axial velocity is: 

a e U e = ^2 a nb U nb + ^ + ( P P~ P &> A e (2.29) 

where the change in pressure with axial distance has been 
pulled out of the source terms. A e is the area term on 
which the pressure is acting. A shorthand expression is 
used for neighbor coefficients and velocities. A very 
similar discretization equation for radial velocity can be 


formed. 
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In the pressure solution process, the solutions of 
velocity and pressure is divided into two components. The 
current velocity and pressure solutions are denoted if, V 
and P* . Due to an inexact solution there is a pressure 
correction P' which causes corresponding velocity 
corrections U' , V' . 


u = u* + u‘, 

(2.30a) 

V = V* + V 1 , 

(2.30b) 

p = p m + p' 

(2.30c) 


A velocity correction equation may be formed by 
plugging in if and V* into the discretized momentum 
equations. Subtracting the equation for U* from the 
equation for U gives a velocity correction equation: 


< 2 - 311 

The first term is unknown, but, disappears at convergence, 
so is ignored. This gives: 



(2.32) 


Plugging this in the equation for U gives the velocity 
correction equation: 

a. = dJ ♦ ^ (2.33) 


The continuity equation is integrated over the control 
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volume for node P. The resulting equation involves 
velocities at the faces of the control volume. The 
equations for velocity correction are then substituted into 
this equation. A discretized for P' equation is formed by 
rearranging the terms : 

apPp = a^p's + atp'„ + A,p' N + a^s + b, ( 2 . 34 ) 

where b includes the mass imbalance due to if and V* and the 

change in density with time. 

This procedure for the solution of convection-diffusion 
problems is a very brief description of the Semi- Implicit 
Method for Pressure Linked Equations. The SIMPLE algorithm 
and its variants are the most used solution methods for low 
Mach number combustor flows. A much more general and 
thorough development of the SIMPLE method is given in 
Patankar [1] . 

Computational models using the SIMPLE algorithm and 
hybrid differencing have been widely used and have proven 
robust for combustor flows. The accuracy of the hybrid 
differencing scheme is believed to deteriorate when the flow 
does not align with the grid. Correa and Shyy et al. 8 have 
studied second order upwind differencing and the QUICK 
differencing scheme of Leonard 9 . Using simple simulations, 
it was demonstrated that one scheme was not generally 
superior for all cases. Higher order schemes are generally 
thought to have problems with robustness. However, Correa 
and Shyy stated that convergence rates of the more accurate 
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numerical schemes were comparable to the hybrid scheme in 
the context of curvilinear coordinates. It was also found 
that with an appropriate grid distribution, the difference 
between finite difference operators could be quite small. 

2.4 - Grid System and Axisymmetric Equations 

The grid system used in this work is orthogonal. The 
grid setup is Practice B cited by Patankar [1] . Practice B 
first sets up the control volumes and then places the scalar 
gridpoints in the center of the control volume. The method 
used in this work uses plain averages rather than the more 
exact geometric interpolations or harmonic averages for 
physical quantities such as viscosity, etc. at cell faces. 
The boundary control volumes are of infinitesimal thickness, 
which eliminates half control volumes for scalar values. 
Axisymmetric geometries are solved in this thesis. The 
equations are slightly different than for Cartesian 
geometries. The general equation for orthogonal 
axisymmetric geometries is: 


a 

at 


<p4>) + 



(2.35) 


The source terms and diffusion coefficients for the various 
momentum and scalar equations are given in Appendix 1. 

2.5 - Boundary conditions 

Inlet conditions are read in and set the main program. 
Initial values are set in subroutine INIT. A restart 


lil IT 
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capability is also in the main program. Each variable has 
its own calculation subroutine and entry point in the 
boundary condition subroutine. The Boundary conditions are 
implemented in subroutine MODPRO. The program must be 
modified for different geometries. 

2.5.1 - Inlet conditions 

All inlet conditions are specified. A large portion of 
this work is to validate computational algorithms. This 
requires that inlet conditions be as well specified as 
possible, so that quality of numerical predictions can be 
judged. Inlet velocities, temperatures, species mass 
fractions, and turbulence kinetic energies are specified 
based on actual measurements, if possible. If there are 
multiple inlets, the inlet pressure is set to be the same 

for both. 

2.5.2 - Axis of Symmetry 

Radial velocity is zero at the axis of symmetry. The 
radial gradient of all variables, except radial velocity is 
zero at the axis of symmetry. Mathematically: 

= 0 

\ A^/r=0 

The axis of symmetry is located on the j=l gridline, 
axisymmetric boundary condition produces the result: 

<j> . J= i = except for velocity V 

This boundary contribution is implemented each iteration. 
During the iteration process the neighbor coefficient for 


(2.36) 
The 

(2.37) 



the axis of symmetry is zero as the area term is zero. 

2.5.3 - Outlet 

At the exit of the combustor the axial gradient of all 
variables except for axial velocity is assumed to vanish. 
Thus, the neighbor coefficient a E is set to zero at the 
exit. Thus, exit variables, except axial velocity, need not 
be specified. 

For axial velocity, an exit velocity correction due to 
mass conservation is applied every iteration. Generally, 
this greatly helps the rate of convergence. 

2.5.4 - Solid walls 

The law of the wall is used for velocities parallel to 
solid walls. The laminar sublayer extends to y + of 11.63. 
The value of y + is calculated by using the density, laminar 
viscosity, turbulent kinetic energy and appropriate length 
for the grid cell adjacent to the wall. Thus for axial 
velocity: 


pC^VE Ay/2 ^ 


(2.38) 


If y* is less than 11.63, the shear stress is laminar and is 
given by: 



(2.39) 


Otherwise the cell is in the turbulent flow regime and shear 

stress is calculated using the law of the wall: 

where K is taken as 0.4187. Shear stresses are included in 
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In (9y + ) 


(2.40) 


the momentum equations as source terms. 

For turbulent kinetic energy, the neighbor coefficient 
corresponding to the wall is set to zero and the source term 
is : 

S* = \i c G - [2c£ /4 pic 3/2 y V Ay], y + < 11-63 (2.41a) 

S k = - [2C,? /4 pJe 3/4 ln (9y*) / (KAy)], y + ^11.63 (2.41b) 

The turbulent dissipation rate next to the wall is 
calculated by assuming that the rates of generation of 
turbulent kinetic energy and dissipation are equal so that 
the value of £ is: 


e 


2c£ /4 k 3 /2 

KAy 


next to the wall 


(2.42) 


The value of £ next to walls is fixed by the forcing the 
source term coefficients to overwhelm all other neighbor 
contributions : 


S p = -1 . 0 x 10 30 , 


(2.43a) 


S c = 1 . 0 x 10 30 e 


(2.43b) 


Walls are treated as adiabatic. Thus, the neighbor 
coefficients for wall influences for heat transfer and 
diffusion of species is set to zero. 

2.6 - Spray Combustion Modeling 
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Spray modeling has been under development for about 
three decades. Many review type papers on spray combustion 
modeling have been written. These include, Faeth' s paper in 
1977 10 , and Law's paper in 1982 11 , which cites improvements 
to original droplet evaporation modeling. The 
sophistication of spray models being used for combustor 
calculations is increasing. The earliest spray models 
treated spray vaporization as a single droplet in quiescent 
flow. The difference in velocity between the droplet or 
fluid phase and the gas phase velocity was ignored. In a 
1983 12 paper, Faeth cites models which include 
consideration of the difference between the droplet velocity 
and gas phase velocity. In Sirignano's paper in 1983 13 , 
fluid flow within the droplet itself was incorporated in a 
spray model. 

One of the earliest studies of spray combustion was 
performed by Onuma and Ogasawara 14 . Some conclusions were 
that for the spray burner used, the region where the 
droplets exist is limited to a small area above the burner 
nozzle, and that it was reasonable that droplets do not burn 
individually, but, that the vapor cloud formed by the 
vaporization of the droplets burns like a diffusion flame. 
Many past calculations have successfully ignored spray 
vaporization or the fluid phase. Fuel was treated as being 
completely in the gaseous phase. However, as the demand for 
more efficient and lower emissions combustors has increased 
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along with the needed computer resources, spray modeling is 
being developed and implemented in the prediction of gas 
turbine combustors. 

According to Law, the basic droplet combustion model 
was formulated by various people, including Spaulding 15 . 

This model has since been termed the d *- Law because the 
model predicts a linear relationship between the square of 
droplet diameter and time. Many improvements have been made 
to the cP-Law. These include modeling of droplet heating, 
variable specific heat in the gas phase, and the multi- 
component fuels. 

Actual sprays consists of a fairly continuous 
distribution of droplet size and droplet number density. 
Spray modeling employs discrete droplets in the prediction 
of spray combustion. Onuma and Ogasawara [14] thought that 
modeling a single characteristic droplet would be sufficient 
to accurately calculate spray combustion. The 
characteristic droplet is determined by some sort of 
averaging process regarding initial droplet size, position 
and velocity. The effect of multiple droplets or groups is 
found by taking multiples of the characteristic droplet 
effects. Most gas turbine spray modeling assumes that 
droplets do not interact, that the spray is dilute. Thus, 
the calculation of spray vaporization is treated as a 
function of individual droplet characteristics and the 
surrounding gas phase. Current models track numerous 
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droplet groups. A range of different size, position and 
velocities may be used to approximate a droplet 
distribution. The most computationally expensive droplet 
models even model turbulent droplet dispersion. This is 
done by incorporating the effect that a random turbulent 
eddy would have on the droplet. Enough droplets must be 
computed to accurately predict overall spray. 

Consideration is being given to modeling non-dilute or 
dense sprays, which exist in fuel atomizers. Faeth's paper 
in 1987 16 includes a section on the modeling of droplet to 
droplet interactions or dense sprays. Tsai and Sterling 17 , 
show that the interaction between droplets increases as the 
distance between droplets decreases. There is difficulty in 
measuring dense spray properties. Non- intrusive optical 
methods are usually used. When multiple droplets are 
detected in a measurement volume, the measurement is 
typically invalidated. 

The fuel break up process is non-uniform. In a section 
on the design and development of gas turbine combustors, 

A.H. Lefebre 16 writes that the fuel atomization is normally 
accomplished by spreading the fuel into a thin sheet and 
then increasing its surface area until it becomes unstable 
and disintegrates into threads or ligaments. These threads 
then break up into rows of droplets. The ligaments and 
droplets are highly nonuniform. The droplets are of various 
sizes and shapes and are not generally spherical. Droplet 
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vaporization modeling assumes perfectly spherical droplets. 
It simply is too difficult to calculate non-spherical 
droplets in combustors or measure them for that matter. The 
break up of fuel stream into droplets is generally handled 
by empirical functions, rather than any physical modeling. 
2.7 - Spray Solution Method 

A spray combustion model including most of the above 
mentioned characteristics was modified. Simplified internal 
droplet fluid flow and temperature distribution is handled 
by the simplified vortex model of Tong and Sirignano 19 . 

The spray model is a simplification of the multicomponent 
spray modeling done by Raju and Sirignano 2 . The model 
assumes a dilute spray. Constant properties are assumed 
within the droplets. The droplet trajectories are solved 
using discretized Lagrangian equations. The turbulent 
dispersion of droplets isn't considered, as it would add at 
least an order of magnitude time to spray computations. 

Sprays consist of many thousands of droplet particles 
each having their own velocity, drag forces, position, 
droplet diameter, and distributions of temperature and 
vortex strength. Due to limited numerical capabilities, a 
limited number of droplets are calculated. Each 
characteristic droplet has its own axial position x k , and 
axial velocity U k , which is a function of time. The 
subscript k refers to the kth droplet group. The change in 
droplet position with time is given by: 
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dx k 

~dt " Uk ' 


(2.44) 


= V 

dt V *' 


dz k 

dt 




(2.45) 


(2.46) 


The droplet velocities are affected by drag effects caused 
by differences in velocity with the surrounding gas phase: 


dt 


dt 
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(2.47) 


(2.48) 


(2.49) 


Where the subscript g denotes the gas phase, r k denotes the 
radius of the droplet and Re k is the droplet Reynolds 
number taken with respect to the difference in velocity with 
the surrounding gas phase: 


Re k = 2 OT g -U k ) 2 + (V g -V k ) 2 + ( W g -W k ) 2 ] 1/2 (2.50) 

P s 


The coefficient of drag is a function of droplet Reynolds 
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number: 
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_ 24 f . tekl) (2.51) 

D Re k { 6 ) 

Molecular viscosity is evaluated by Sutherlands ' s equation: 

M” - 1.4637x10^-5^5 < 2 - 52) 

A reference temperature is used for calculation of viscosity 
in the vicinity of the droplet: 


Tief 


T r - + 



(2.53) 


where the subscript . s refers to the droplet surface. 

The droplet regression rate is a composite function of 
droplet Reynolds number, Re k . Three ranges are used to 
approximate the function: 


d(r 2 ) 

1 

3.| 

CM 

1 

II 

—Re k 

U2 f (B k ) , 

Re k > 20 

(2.54a) 

dt 

Pjcl 

71 





d(r±) M iLlfi + ( i +Re k ) In ( 1 +B k ) , KRe k <l 

dt p* L 


(2.54b) 


d tip- = it [i + ( i +Re k ) 1/3 ]ln ( l +B k ) , Re k < l 


where B k is the Spaulding transfer number. The function 


f(B k ) is obtained from the solution of Emmon's problem. 


The internal droplet temperature is calculated from: 


dT k 

dt 


17- 


c piPi r k L a « a 


, , . . dT k 

a-^-£ + (1+C< t) a) * 


da 


(2.55) 


where Jcj refers to thermal conductivity in the droplet or 
liquid phase, and C(t) is given by: 


C(t) 


3 ^jiPi 
17 iCj 



(2.56) 


a represents the streamline of a Hill's Vortex in the 
droplet. The boundary conditions for the droplet 
temperature are: 


= _L 

da 17 



a=0 (2.57a) 
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7 d(r 2 ) 

±k dt~ 


a = 1 (2.57b) 


where l k is the latent heat of vaporization of the fuel. 
The initial droplet temperature is set to the droplet 
injection temperature. The Spaulding transfer number is: 


Bk . 


L Jc, eff 


(Yfs-Yf) 

<1 -Yfs) 


(2.58) 


Ik.eff is the effective latent heat of vaporization modified 
for heat loss to the droplet interior: 
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Jr, err 



(2.59) 


where m k is the droplet vaporization rate. The fuel mass 
fraction at the droplet surface is given by: 

r,. - < 2 - 60 

where W a is the molecular weight of the gas excluding the 
fuel and x is the m °l e fraction which is found by assuming 
phase equilibrium at the droplet surface, and using the 
Clausius-Clapeyron relation: 



where P n is the normal atmospheric pressure, P is the 
prevailing pressure, l a is a constant and T b is the boiling 
temperature for the droplet. More accurate expressions than 
the Clausius-Clapeyron relation involving additional 
parameters can be found in Reid et al . 

The spray provides a mass source term in the continuity 
equation. The droplets also affect the momentum equations 
through drag forces. Thus, appropriate source terms must be 
included in the gas-phase equations. The source terms are 
listed in Appendix 2. 
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Chapter 3 

NUMERICAL SIMULATION OF COMBUSTION 
The field of combustion demands knowledge in compu- 
tational fluid dynamics, turbulence modeling, chemical 
kinetics modeling of combustion, real gas effects, combustor 
concepts and turbulent-chemistry effects. A variety of 
combustion models exist. Combustion models cannot be 
implemented as easily as changing damping coefficients for a 
turbulence model. Many of these turbulent combustion models 
are application specific, that is, the combustion model was 
developed for a specific combustion regime. The combustion 
regime is decided by the length scale of combustion, which 
is dependent on the rates of mixing and the rates of 
chemical reactions. Much combustion model development of 
late has involved combustion for jet flames where the com- 
bustion scale is of the order of the smallest scales of 
turbulence. The chemical reaction rate is very fast 
compared to the mixing rate. Unf ortunately, combustion is 
distributed in practical combustors due to the use of slower 
reacting, more complex fuels and high turbulence levels. 

The expertise in computational combustion is divided 
into two main groups. One group very closely models 
thermodynamic and transport properties. Turbulence is 
generally ignored by this group. If there is transport, it 
is laminar, and species diffusive velocities, viscosities, 
etc. are precisely modeled. Equilibrium, partial 
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equilibrium, or complete system of chemical reactions may be 
solved. Equilibrium assumes infinite time for chemical 
reaction. Species concentrations do not change after 
reaching equilibrium. Reactions are still proceeding, but, 
forward reactions are perfectly balanced by reverse 
reactions. The final chemical species for equilibrium 
combustion calculations are found by minimizing the 
Helmholtz or Gibb's free energy depending on whether it is a 
constant-volume or constant -pressure process. Finite rate 
combustion involves time stepping through a large number of 
ordinary differential equations, which are usually stiff. 
Small time steps must be used to retain numerical stability. 
The combustion process is viewed as a collection of a number 
of elementary reactions involving likely combinations of 
reacting species. Each of these elementary reaction steps 
has its own constants for activation energy, pre-exponential 
constant, the species concentration dependence and possibly 
some extra temperature dependence. Many of these elementary 
reactions are reversible and many of the species are in 
numerous competing chemical reactions. Some reaction rate 
constants may be calibrated using shock tube experiments. 
Other reactions are currently unknown, so combustion rates 
are inferred. When the kinetic scheme duplicates shock-tube 
data over a range of conditions the scheme is said to be 
accurate or validated. Not too surprisingly, there are 
multiple schemes which have been semi— validated for the 



38 


combustion of simple fuels. Improving the accuracy of the 
model is done by adding more reactions, which increases the 
applicability of the kinetic combustion model over wider 
temperature, pressure and species concentration ranges. 

There are multiple competing paths in which fuel and oxidant 
are burned. The influence of these competing reaction paths 
is highly affected by temperature. At different temperature 
ranges, different reaction paths can take precedence. 
Including multiple reaction paths extends the validity of 
the kinetic model. The current state of chemical kinetic 
modeling is that the combustion of hydrogen and some of the 
lower hydrocarbons is quite accurate. Laminar flame speed 
and heat release are well predicted for the very simple 
fuels. The kinetic combustion modeling of heavier fuels 
used in practical combustors, such as kerosene, is still in 
a development stage. Kinetic modelers will not advocate 
combustion schemes for heavier fuels until the combustion 
models for the lighter fuels are sufficiently validated, as 
the heavier fuels break down into lighter fuels and other 
intermediate species during the combustion process. 

The fine degree of kinetics combustion modeling does 
not lend itself to calculation of multi-dimensional 
turbulent combustion flow. They are simply too 
computationally expensive. Thus, the implementation of 
kinetic reaction submodels into multi-dimensional CFD models 
has only involved hydrogen or methane with no turbulent- 
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chemistry interaction. That is, mean quantities are used in 
the reaction rates, there is no turbulence correction. 
Currently, the solution of species transport for 
recirculating flows using CFD is such that the calculation 
of individual species without chemical reaction is 
inaccurate. That is, the transport of one species can 
easily be off by over ten percent, a level of accuracy CFD 
calculations should meet to have predictive validity in the 
design of low-emission combustors. Designers would like to 
be able to predict within a few percent. Various combustion 
models have been implemented to predict combustor flows in a 
reasonable time. 

3.1 - Infinite Rate Combustion Models 

The simplest CFD combustion model is the fast or 
infinite rate reaction model. This type of scheme does a 
good job of predicting the exit temperature profiles of 
combustors. In this scheme, once the fuel and oxidant are 
mixed, the reaction rate is such that the reactants are 
immediately converted to products. The reactants are 
consumed in strict proportions. A simple example is the 
stoichiometric burning of methane. The reaction is referred 
to as stoichiometric as it assumes stoichiometric burning of 

CH a + 2 0 2 - C0 2 + 2 H 2 0 

the fuel and oxidant. One molecule of methane is consumed 
with two oxygen molecules to produce one molecule of carbon 
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dioxide and two of water. In units of mass 16 Kg of fuel 
react with 64 Kg of oxygen. Four Kg of oxygen react with 
each Kg of methane. The reaction rate source terms are 
related. 

Rfu = R ox /* ( 3 - 1) 

Using the general transport equation for species, 
equations may be written for the mass fractions of oxygen 
and fuel. Subtracting a constant specific portion of the 
transport equation of oxidant to that of fuel eliminates the 
chemical reaction source term, leaving a conserved scalar. 
Using a conserved scalar eliminates calculating troublesome 
reaction rates. In the previous example, if we form a 
variable composed of the fuel mass fraction minus one 
quarter of mass fraction of oxygen the troublesome reaction 
rate term is eliminated. The general equation for this 
conserved scalar is: 

f = Y f -(F/0) st Y 0 (3.2) 

The term (F/0) st is the fuel to oxidant molecular weight 
ratio for stoichiometric burning process. Other conserved 
variables could be formed. Another popular conserved 
variable is the sum of all fuel fractions, both unreacted 
and reacted. This conserved variable is often referred to 
as the mixture fraction. The solution of the conserved 
variable defined in (3.2) gives the fuel and oxidant mass 
fractions. If the conserved variable is positive, the fuel 
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mass fraction is the same as the conserved variable, and the 
mass fraction of oxygen is zero. If the conserved variable 
is negative, the fuel mass fraction is zero and the mass 
fraction of oxygen is 4 times the negative of the mixture 
fraction. The sum of the mass fractions of all species is 
one so the mass fraction of products is simply. 

Y pr = 1.0 ~ Yy - Y 0 (3.3) 

An atom balance is done to figure the individual product 
species fractions. If the oxidant is air, the above 
equations must be slightly modified to include nitrogen. In 
this infinite rate combustion model, the width of the 
reaction zone is zero. Combustion occurs in a flame sheet. 
There is no overlap of oxidant and fuel. Fuel exists on one 
side of the combustion sheet, and oxygen on the other. 

An enthalpy equation is solved from which temperature 
can be found from the various mass fractions. Radiation is 
typically ignored as research combustors are usually 
operated at lower temperatures and pressures. The equation 
of state is used to solve for density: 

p =PM/R U T ( 3 • 4 ) 

The equation uses the universal gas constant and the mixture 
molecular weight (inverse of sum of number of moles of each 
species divided by each species molecular weight) . The 
mixture molecular weight in terms of mass fractions is . 
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Various improvements can be included, such as allowing for 
variable specific heat with temperature. As specific heat 
increases with temperature, this change alone can improve 
the prediction of temperature by a couple of hundred 
degrees. The use of infinite reaction rates is known to 
give good predictions of exit combustor temperature 
profiles, and major species concentrations. Partially 
reacted species or nitrous oxides are not predicted by 
infinite rate combustion models. Rather than assume 
complete chemical reaction, equilibrium combustion can be 
modeled. At equilibrium the species concentrations are 
fixed, but, the rate of forward chemical reaction is just 
balanced by the rate of back or reverse chemical reaction. 
Equilibrium modeling allows for the presence of minor 
chemical species, allows reactant species to overlap and a 
broadening of the reaction zone width. Carbon monoxide is 
usually overpredicted with rich fuel mixtures. Equilibrium 
combustion models significantly overpredict chemical 
reaction and product species at low temperature. Formation 
of some species such as nitrous oxide is a function of time. 
The equilibrium model does a very poor job of predicting 
nitrous oxides. 

One model that allows for overlap of reactant species 
at low temperature is the finite rate combustion model. 
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3.2 - Finite Rate Combustion Models 

Finite rate combustion means that fuel and oxidant can 
co-existing at the same time. A conserved scalar is still 
conserved in finite rate combustion modeling, but, the 
concentration of fuel or oxidant cannot be solely obtained 
from the conserved variable as before. An additional 
equation must be solved for that includes a finite rate 
reaction source term. The calculation of the finite rate 
reaction term is usually based on arrhenius reaction rates. 
Arrhenius was the first to include the Boltzmann factor into 
chemical reaction rates. At room temperature oxygen and 
hydrogen can co-exist with no reaction. Energy must be 
supplied to begin the combustion process. Arrhenius 
postulated that collisions between fuel and oxidant do not 
automatically produce combustion. The collision between 
molecules must be strong enough that an energy barrier is 
overcome. The Boltzmann factor takes this into account. 

The equation: 

k = exp(--^,) (3 - 6) 

is called Arrhenius law. A typical Arrhenius reaction rate 
may look like 

-^H=-A(C fu ) a (C ox ) i) exp(-^) (3 * 7) 

where C is the molar concentration. By employing the 
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perfect gas law, the expression may be algebraically 
manipulated to use mass fractions which is more commonly 
used in CFD. A is the frequency number to take into account 
the collisions between molecules. Sometimes a dependency on 
temperature to some power is also included. For fundamental 
or elementary chemical reactions, the exponents a and b are 
whole numbers from the molar coefficients in the chemical 
reaction equation. Fundamental chemical equations are 
reactions that actually occur. The combustion of hydrogen 
with oxygen is. globally represented by: 

H 2 + 1/2 0 2 - H 2 0 

The actual combustion process occurs as a series of 
reactions of which only a few are: 

OH + H 2 - H 2 0 + H 
H + O z -* OH + O 
O +H 2 -* OH + H 
OH + OH + M - O + H 2 0 
OH + H0 2 - 0 2 + H 2 0 
H + 0 2 + M - Af + H0 2 

There are numerous elementary reactions to get the final 
product, H 2 0 . In the process there are many reaction 
radicals which are also in many elementary competing 
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reactions. To further complicate things, the intermediate 
species, H 2 0 2 , is also important. R.J. Kumar 22 , and S. W. 

Kim 23 have successfully used a set of 24 elementary 
reactions involving 9 species to predict hydrogen- oxygen 
combustion. This set of elementary reactions was reduced 
from a set of 48. At different temperature ranges different 
reaction paths take dominance. There is wide disparity m 
the reaction rates which causes great difficulty in their 
solution. Some kinetic reaction schemes are validated by 
duplicating calculated flame structure, laminar flame 
speeds, and chemical composition using large numbers of 
fundamental reactions. Computational models used to do this 
are one - dimens i ona 1 , include modeling of individual species 
laminar diffusion velocities and use very small time steps. 

In multi-dimensional recirculating flow global reaction 
models are typically used. The global reaction models are 
only calibrated over limited temperature, species 
•concentration and pressure ranges. In a 1972 paper, Dryer 
and Classman 24 used a temperature range of 1030 to 12 30 K 
for oxidation of carbon monoxide and methane. The most 
often cited global reaction schemes are from Westbrook and 
Dryer 25 . The reaction schemes were designed to capture 
rich and lean flammability limits, information on flame 
temperature and burned gas composition and flame speed. 

They do not accurately describe the chemical structure of 
Different models can produce different results 


the flame. 
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with the same reaction numbers. Coffee 26 recommends higher 
frequency numbers as his results were slightly different 
than Westbrook and Dryer's, unless he adjusted the reaction 
numbers. There is further controversy in regards as to what 
phenomena should have primary importance in modeling. 
Coffee 27 , in another report, states that heat release 
should be used to calibrate reaction numbers, not the 
criteria used by Westbrook and Dryer. The frequency number 
given in this paper is orders of magnitude higher. 

If one wants to calculate intermediate species such as 
carbon monoxide, at least one additional equation must be 
solved. Nikjooy et al. 28 gives an example with the 
variables and source terms . A four step scheme has been 
used by Srivatsa 29 . This scheme is from Hautman et al. 30 
and allows for additional species of Hydrogen and ethylene. 
This scheme includes product concentration in the reaction 
rates. The scheme was validated over 960 to 1540 K for 
propane . 

An early study of global reaction schemes was done by 
Abdalla et al. 31 One conclusion was that "inaccuracies 
could arise from the unrealistic simplification of chemical 
kinetics and the variable effects of turbulence" . Duterque 
et al. 32 also came up with a global scheme for methane, 
propane, benzene and isoctane validated for a well-stirred 
reactor model. The authors stated that the model could be 
used in the predictions of burners and combustors. 
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In a more recent paper, Westbrook and Dryer 33 discuss 
detailed kinetic modeling for various hydrocarbon fuels 
along with a short section on global modeling. They note 
that using a single reaction scheme leads to overprediction 
of products and temperature. In 1986, Paczko et al. 34 gave 
a four-step reaction mechanism for methane and a five step 
reaction mechanism for propane. The propane mechanism 
includes two reactions for the initial break up of propane. 

A three-step mechanism is recommended by Seshadri and 
Peters 35 for stoichiometric methane-air flames. The 
mechanism was tested for a counterflow diffusion flame. In 
1987, Kiehne et al. 36 discuss an improved eight-step 
mechanism for propane. The model is based on the original 
work of Hautman et al. The model better predicts ignition 
delay and combustion at higher pressures. In 1988, Jones 
and Lindstedt 37 proposed a four-step global mechanism for 
the combustion of alkane hydrocarbons up to butane. The 
scheme was tested on methane and propane flames, along with 
diffusion flame data for a methane-air flame. Good 
agreement was predicted for flame speed, flame thickness, 
and species profiles. Bilger and Starner 38 constructed a 
four-step mechanism from a set of 58 elementary reactions. 
The model uses steady- state approximation for the oxygen 
atom rather than a partial-equilibrium approximation. The 
mechanism was used to predict flame structure for the 
counterflow diffusion flame. Good predictions were obtained 


for minor species over a range of flame stretch. 

Some of these schemes involve the implementation of 
source terms for calculating product species. This 
researcher found problems when calculating product species 
for multi -dimensional recirculating flows. The frequency 
numbers used in arrhenius rates are fairly large. At higher 
temperatures, arrhenius reaction rates can predict rapid 
reaction even though reactant concentrations are quite low. 
The arrhenius reaction rate is only accurate for an 
extremely short time step, not the typical time steps used 
in implicit calculations. Thus, product concentrations can 
be overpredicted at high temperatures. Most often, 
overpredicted concentrations are simply overwritten during 
the iteration process. 

3.3 - Turbulent Combustion 

All practical gas turbine combustors operate in the 
turbulent regime. Thorough mixing and burning must occur 
within rather short axial distances in current gas turbine 
combustors. Turbulence for chemical species means that 
there is variation with time in the chemical species and 
temperature. A consequence of this is that even with 
infinite reaction rates, fuel and oxidant can exist at the 
same location, but, not at the same time. Where this 
occurs, the time averaged concentration of fuel and oxidant 
will not be zero. Many different turbulent combustion 
models have been worked on in the last two and a half 
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decades. Many of these models do well within the range of 
their assumptions. 

3.3.1 - Eddy Break Up Combustion Models 

One of the earliest turbulent combustion models is the 
eddy-breakup model proposed by Spaulding 39 . The original 
application for the model was for premixed combustion. It 
was observed that changes in velocity and fuel to air ratio 
had little affect on combustion in premixed combustors. The 
use of so-called laminar arrhenius reaction rate combustion 
models didn't predict this. The controlling mechanism 
wasn't temperature related, but, was controlled by mixing 
rates. In the eddy-breakup model combustion is viewed as a 
mixture of burned and unburned eddies. The combustion 
process is controlled by the mixing and exchange of energy 
from burned eddies to neighboring unburned eddies. The 
eddies are stretched and folded around each other until 
local gradients are strong enough that mixing occurs. This 
is very similar to the cascade of turbulent kinetic energy 
from large eddies to smallest eddies where the energy is 
finally dissipated. It was further proposed that combustion 
was also proportional to reactant concentration. If oxygen 
concentration was low, as if most of the oxygen were already 
consumed, the reaction rate should be lower. The initial 
combustion model did a good job of qualitatively predicting 
the correct hydrodynamic phenomena for premixed combustors . 
The eddy-breakup model was modified over the next few years. 
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In another paper 40 , it was proposed that reaction rates 
should be proportional to the concentration fluctuation not 
the reactant concentration. An equation for concentration 
fluctuation was proposed and successfully tested for a 
turbulent jet. In another paper. Mason and Spaulding 41 
tested the combustion model. From a later review 42 , it was 
reported that the improved more physically realistic model 
predictions were worse than previous modeling. In 1976, 
Magnussen and Hjertager 43 proposed an another eddy-breakup 
model where the reaction term was proportional to the 
minimum of related constants times appropriate fuel, 
oxidant, or product concentration. The model performed well 
results, but, the constants had to be increased for a 
different case. Isothermal turbulence modeling was also 
developing. The k-e turbulence model was developed and has 
proven successful. Using the k-e turbulence model, the form 
of Spaulding's eddy-breakup model is: 

~S P = (3.8) 

A measure of the rate of eddy breakdown is the reciprocal of 
eddy lifetime. A measure of eddy-lifetime is the ratio of 
turbulent energy, k, over the dissipation rate, e. 

The eddy-breakup rate was developed for premixed 
combustors but has been applied to other combustors. In 
work at Garrett research 44 the eddy-breakup model was 
combined with an Arrhenius type reaction rate. The reaction 
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rate was taken as the minimum of the Arrhenius rate or the 
eddy-breakup rate. Multiple step chemical reactions were 
also implemented. The combination of eddy-breakup model and 
arrhenius model is easy to implement. The model was 
superior to earlier models. 

3.3.2 - Probability Density Function Models 

The eddy-breakup model has been long regarded as a 
semi-empirical combustion model, even though it has proven 
successful. In turbulence modeling, the current thought is 
to include as much of the physical process as possible to 
improve the models. Thus, there has been great effort to 
model various turbulent terms in the averaged Naviar Stokes 
equations for chemical species. One method is to 
incorporate a correction factor for turbulence which can be 
applied to the laminar Arrhenius reaction rate. Employing 
Reynolds decomposition and the perfect gas law, a typical 
turbulent reaction rate may look like: 


w fu~ 


-BP 2 


(Y f +Y'„) (Y c 


- 7 ? 


(3.9) 


Using T a = EJR^ and some algebraic manipulation, the 
exponential part of the above expression is: 
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The second exponential may be expanded in an infinite 
series. This infinite series is only convergent if T a is 
not large compared to the averaged temperature and if the 
variation in temperature is not large compared to the 
averaged temperature. This is not probable in practical 
combustors. An additional problem with this method is 
modeling the correlation term Y' fu Y' ox . One way around the 
reaction rate closure problem is to treat turbulent reacting 
flow as a stochastic process 45 . That is, the solution of 
average scalar quantities involves the weighting of the 
scalar with the probability of the scalar. 

1 = f 1 p(f) df ( 3 . 11 ) 

Jo 

77 = [ 1 Y i (f)P(f)df ( 3 . 12 ) 

Jo 


Here it is implicitly understood that the probability 
density function is also a function of time and position. 
The first stochastic methods used in multi-dimensional 
turbulent combustor flow involved the use of two-parameter 
probability density functions. The pdf is completely 
specified by two parameters, the mean and variance of a 
conserved scalar. These are referred to as the first and 
second moments. 

One of the basic assumptions of this model is that the 
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reaction is very fast or near equilibrium. This allows the 
use of a single variable, mixture fraction, to specify the 
mixture. Physically, the fast reaction assumption means 
there is little or no overlap of reactant species at any 
instance in time. The incorporation of stochastic modeling 
of turbulence allows the overlap of reactant species over 
time. Much like the eddy-breakup model, reactants occur in 
separate parcels. Over time these parcels move around so 
that the time average of reactant concentrations overlap. 

This occurs with mixtures near stoichiometric conditions 
causing a broadening of the reaction zone. As actual gas 
turbine combustion zones are distributed, this modeling 
better portrays actual combustion. 

Typical two parameter probability density functions are 
rectangular or battlement pdf, Gaussian pdf, clipped 
Gaussian pdf and Beta pdf. Initially, rectangular or 
battlement probability density functions were used due to 
their simplicity. Gaussian pdf's are more appropriate in 
highly turbulent flows, but, are more complicated to 
calculate. Gaussian pdf's are determined by their mode or 
most probable value and standard deviation, not the mean and 
variance which is found by solving their CFD transport 
equations. There is no direct way of solving for the 
Gaussian pdf from the scalar mean and variance. The pdf 
must be specified so that other scalar averages, such as 
temperature and density can be solved. The most efficient 
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way of solving for the Gaussian pdf is to construct a table 
lookup of mode and standard deviation versus mean and 
variance. Thus, when the mean and variance of conserved 
scalar are found during the iterative process, an 
interpolation of the table is done to find mode and standard 
deviation which specifies the gaussian pdf. 

One problem with Gaussian pdf's is the tails which 
extend to infinity. These tails can correspond to 
physically impossible states. The clipped gaussian was 
developed to take care of this. Delta functions are added 
to take the place of clipped tails. The solution of clipped 
gaussian pdf's is slightly more difficult. A pdf which 
resembles the gaussian pdf is the Beta pdf, which was 
proposed by Richardson 46 . The Beta pdf has been 
successfully used by a number of researchers 47 . Unlike the 
Gaussian pdf, the Beta pdf can be directly calculated from 
the mean and variance. 

The conserved variable or mixture fraction is typically 
non-dimensionalized for calculations. The mixture fractions 
at the fuel and oxidant inlet streams represent the range of 
mixture fraction in the combustor. 

o [ Y F ~ (F/0) sC Y o ] m - \Y f - (F/O) st Y 0 ] a ^2 13 ) 

° = [Y f - ( F/0) sC Y o ] f - [Y f - ( F/0) st Y Q ] A 


The mean mixture fraction is found by solving the averaged 
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equation for the mixture fraction. The variance of mixture 
fraction is often referred to as g: 


g = e»2 = J'(0-d) 2 Pr(d)dO (3 - 14) 

0 

where the overbar refers to Favre or density averaging 
rather than the normal time averaging. A modeled transport 
equation is solved for the local variance. The formula for 

beta pdf is: 


Fix) = 


r ( a+b) 
T(a)T(b) 


( l - x )*" 1 


where Osxsl, 
and a>0, Jb>0 


( 3 . 15 ) 


the formulas to find a and b from mean and variance are 


a 




1 ] , 


( 3 . 16 ) 


The major limitation of this method is that computed 
mixtures are close to equilibrium. Equilibrium assumptions 
tend to overpredict CO concentrations. Another choice in 
combustion models is the laminar flamelet model. In 
reviewing experimental work done by Tsuji and Yamoka 48 , 
Bilger 49 noticed that measurements of temperature and 
species concentration, including CO appeared to be functions 
of mixture fraction for laminar diffusion flames despite 
differing shear rates and positions. That is, there exists 
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a state relationship to the mixture fraction. Thus, a table 
of species concentration and temperature versus mixture 
fraction can be constructed from the experimental data. 

This is very useful in computations as the table can be used 
to find thermodynamic quantities during the iteration 
process. This eliminates calculating reaction rates and 
enthalpy. Instead of assuming infinite rate or chemical 
equilibrium, actual experimental data is incorporated into 
numerical predictions. Jeng 50 used this in his thesis to 
perform buoyant free jet combustion calculations using the 
GENMIX 51 computer program. Jeng also modified an 
equilibrium program to check the experimental concentration 
data. The program was modified for rich combustion. Very 
high carbon monoxide concentrations were predicted, 
characteristic of high equivalence ratio equilibrium 
combustion calculations. 

This thermochemical alternative to chemical 
equilibrium was further developed into the flamelet model by 
Liew et al . 52 In this model, reaction rates are very high 
compared to the mixing rate. An instantaneous snapshot of 
the combustion process would show a flame with finite 
thickness or a flame sheet. The interaction of turbulence 
with this flame surface involves movement and wrinkling of 
the flame surface. In a geometrically similar type 
combustor, laminar diffusion flame data is measured. These 
data are then used to provide unique thermochemical 
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relationships solely as a function of mixture fraction. 

These relationships are then coupled with an assumed shape 
pdf to account for turbulence interaction with combustion. 

The actual flame is viewed as an ensemble of laminar 
diffusion flames. This method has produced substantial 
improvements in prediction in fuel-rich regions of an open 
combustor, particularly in CO concentration. 

The laminar diffusion flamelet idea was further 
developed by Liew and Bray 53 for highly stretched 
combustion. It was found that if flames were substantially 
stretched, the state relationship between thermochemical 
values and conserved scalar was affected. As the flame is 
stretched or instantaneous scalar dissipation is 
substantially increased, it was noticed that temperature 
decreased. Also, the mixture fraction at which maximum 
temperature occurs, increases with increasing stretch. The 
region surrounding this maximum temperature is the region of 
■ greatest chemical reaction. Beyond a certain level of flame 
stretch, the imbalance between local heat loss, local active 
species loss, and smaller heat production due to flame 
stretch leads to flame quenching or blowoff. The 
thermochemical or state relationships now involve two 
variables, mixture fraction and flame stretch or average 
instantaneous scalar dissipation. Conceptually, a library 
of flames is measured or calculated as a function of two 
variables. Turbulence effects are still handled by an 
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assumed shape pdf. However, the pdf is now a function of 
two variables, it is a joint pdf. The shape of joint 
probability density functions are much more difficult to 
deduce. The simplest joint pdf is to assume that the two 
variables are completely independent. That is, the joint 
pdf is a product of pdf for mixture fraction and another pdf 
for scalar dissipation. Liew and Bray used a Beta pdf for 
mixture fraction. They assumed that the pdf for scalar 
dissipation was log-normally distributed at high Reynolds 
number. Equations are then formed for the Favre average 
scalar dissipation and other quantities needed to specify 
the pdf for scalar dissipation. This model allows the 
prediction of oxygen penetration into fuel jets where there 
is large strain. The chemistry is laminar and can be 
calculated very accurately. Drake and Blint 54 found that 
finite rate kinetics have an important effect on flame 
structure. The use of partial equilibrium modeling for 
carbon monoxide concentration in low temperature fuel-rich 
flames was shown to be inadequate. The local effects of 
flame stretching correspond to a turbulent ensemble to 
average the flamelet characteristics. 

Flamelet models are based on a one-dimensional 
approximation of the flame. The curvatures of the various 
flamelets are neglected. Coherent vortices which are 
present in many mixing layers are neglected. Borghi 55 , in 
his review of turbulent combustion modeling, states that 
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f lamelet modeling should be confined to cases where the 
chemistry is very fast and the Reynolds number is not too 
large. Peters 56 did a review of the laminar f lamelet 
approach which describes the method more fully. The method 
is an attractive way of treating the prediction of minor 
species in turbulent flow calculations at reasonable 
computer costs. 

Janika and Kollman 57 used a two-variable pdf scheme in 
their studies of turbulent H 2 -air flames. The two pdf 
variables are mixture fraction and reaction completeness. 
This method assumes partial equilibrium for its H 2 /air 
chemistry. The kinetic mechanism selected consists of fast 
shuffle reactions which are considered to be in partial 


equilibrium: 


H + 0 2 ** OH + O ( Rl > 

0 + H 2 ** OH + H ( R2 ) 

H z + OH *» H 2 0 + H ( R3 > 

2 OH - H 2 0 + O ( R4 > 

and slow three-body recombination reactions. 

H + OH + M - H 2 0 + M ( R 5) 

H + 0 + M ~ OH + M (R6) 

H + H + M - H 2 + M ( r7 ) 


The authors stated that the partial equilibrium scheme was 


not good below 1200 K. To reduce the number of variables, 
combined variables are used. This was proposed by Dixon- 
Lewis et al. 58 . The combined variable for hydrogen is: 


-t Mu M„ -> M m 

Ch = C H + — — c OH + — C D + 4-f 
H} - 2 M oh oh M q ° 2 M H H 


(3.18) 


The reaction rate for the combined variable is independent 
of the fast shuffle reactions (R1-R4) . The source terms for 
the hydrogen combined variable are: 

to#., = ~M h ( 2co 5 + 2<o 6 + 2&) 7 ) (3.19) 


The combined variable varies between its unburned and 
equilibrium values. Thus, a nondimens ional reaction 
progress variable is: 





(3.20) 


Transport equations are used to solve for mixture fraction 
and reaction progress variable. The solution of enthalpy, 
mixture fraction, and reaction progress variable, coupled 
with partial equilibrium relations completely specifies the 
instantaneous thermochemical state. Turbulence is handled 
by a combination of assumed shape probability density 
functions. The mixture fraction and combined variable are 
assumed uncorrelated so that the joint pdf is the product of 
the pdf of mixture fraction and another pdf for the reaction 
progress variable. Transport equations for the mean and 



61 


variance of the mixture fraction and progress variable are 
solved. The transport equation for mean reaction progress 
variable must include a reaction term for the recombination 
reactions. A beta pdf can be used for the mixture fraction. 
Three dirac delta functions can be used for the pdf progress 
variable pdf corresponding to unburned, equilibrium and mean 
of the reaction progress variable: 

P(ti) =c 1 8(r|) + c 2 5(ti-t0 + c 3 6(n-l) (3.21) 

The coefficients of the dirac delta terms are determined 
from simple functions of the mean and variance of the 
progress variable. From Janika and Kollman's 59 1982 paper: 



(3.22a) 



(3.22b) 



(3.22c) 


where favre averaging is actually used. 

The scheme was used to calculate the jet diffusion flame of 
Kent and Bilger 60 . The reaction progress variable varied 
between 0.88 and 1.0 across the calculation at x/D of 20. 
The progress values increased further downstream. Janika 
and Kollman state that the more elaborate modeling used in 
the calculation didn't significantly improve prediction of 
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stable components and temperature. However, the modeling 
did considerably help in the prediction of monatomic oxygen 
and reasonable NO concentrations. 

Correa et al. 61 include a couple of CO reactions in 
the above scheme for calculation of CO/H 2 /-air turbulent 
diffusion jet flame. The calculation method successfully 
predicted experimental measurements of superequilibrium 
concentrations of OH. A year later, Correa 62 published 
another paper where a progress variable was added for CO 
concentration. The partial equilibrium equations used in 
the above models are only valid above at least 1200 K. The 
correlation assumptions made about joint probability density 
functions is questionable. Pope and Correa 63 found that 
joint pdf of velocity, conserved scalar and reaction 
progress variable were correlated. The kinetics models for 
higher molecular weight fuels are much more complicated. 
Including additional variables into the joint pdf for 
breakup of the heavier fuels into lighter components is even 
more questionable and untested. Borghi in his 1988 review 
paper generates partial pdf shapes for some situations. At 
best, these shapes only qualitatively resemble Gaussian like 
profiles. The next step in turbulent combustion 
calculations is to actually remove some of the assumptions. 

Some recent papers were written about the current status 
of Combustion modeling. This include Hukam Mongia's 
paper 64 which covers some of the calculations done at the 
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Allison Gas Turbine Division of General Motors Corp. Correa 
and Shyy [8] covers the combustion modeling effort done at 
General Electric ' s Corporate Research and Development 
Center. Much of this effort was to develop their CONCERT 
code. There are also quite a few texts on combustion which 
include: Kuo's 65 text on the principles of combustion, 
Williams' 66 reference on combustion theory, Bartok and 
Sarofim's 67 book on fossil fuel combustion, Chigier's 68 
reference on energy, combustion and environment, and Oran 
' 69 reference on CFD simulation of reactive flow. 
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Chapter 4 

HYBRID COMPOSITION PDF SOLVER 

Many assumptions have been made in developing the 
previous assumed shape pdf models. Implicit in those 
assumptions was that chemistry is very fast. Even 
correcting this by using partial equilibrium solutions for 
jet flames is only valid at higher temperatures and thus, a 
high degree of reactedness. Another assumption was the 
shape of the pdf. The focus of this work will be to further 
develop and test a pdf method which does not make these 
assumptions . 

For this work, the pdf is initially assumed to be a 
function of time, position, temperature or enthalpy, and 
multiple species concentration. This is referred to a 
composition joint pdf. The pdf is called joint because the 
pdf is a function of more than one variable. The hybrid 
composition pdf uses traditional CFD techniques to obtain 
solutions of velocity, pressure, turbulence kinetic energy 
and dissipation. The composition pdf handles the solution 
of species concentration and temperature. The hybrid 
composition pdf solution is thus a combination of 
traditional CFD methods and stochastic modeling. If the pdf 
is assumed to be also a function of velocity, the solution 
method is vastly different from traditional CFD techniques. 
Thus, the composition pdf represents a good intermediate 
step to complete pdf modeling, which is in its infancy. In 
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1989 Anand et al. 70 used the joint velocity pdf method to 
calculate the turbulent flow behind a backward facing step. 
In this calculation only the pressure field was needed from 
another solver. Preliminary work is being done to calculate 
all quantities, including pressure in the pdf solver. In 
1993, Anand, Pope and Mongia presented predictions for a 
swirling flow. The method did not incorporate a radial 
dependence for pressure, although an externally imposed 
pressure gradient could be incorporated into the equations. 
The mean swirl number for the flow was reported to be 0.09. 

As it has been shown that turbulent combustion models 
are rather simplistic and full stochastic methods are not 
developed, the hybrid pdf represents an excellent testbed 
for testing predictions of turbulence effects on combustion. 

Various pdf transport equations can be written. One 
form for the composition pdf is: 

+ pfl tt 3 a P + * • ■ * ^ “ (4.1) 

- 9.( p <»?!♦* <*>■=♦*>« - pgE 3 W<*«l** (x) 

The terms represent the time rate of change/ mean 
convection, chemical reaction, turbulent convection, and a 
mixing term, respectively. The symbol over P denotes Favre 


averaging . 


p = pp/p 


(4.2) 


3 is the scalar dissipation: 

e i.j a (4.3) 

where D is the diffusion coefficient. The use of favre or 
density weighing eliminates some terms in the averaged 
equation. The notation <x| y> denotes the mathematical 
expectation of a random function x conditioned on y. 

Terms on the left side of the pdf evolution equation do 
not need modeling. That is to say, they don't have 
additional unknown variables which demand assumptions or 
relations to produce closure models. The time rate of 
change and mean convection are seen in traditional CFD 
modeling. The last term on the left side is the evolution 
of the pdf in transform space due to chemical reaction. 

This term does not employ modeling assumptions about various 
turbulence quantities such as Y ± ' and T’ . The pdf method 
can use the laminar kinetic reaction rate modeling. As was 
mentioned previously, using regular moment methods to solve 
for chemical reaction produces correlations that have to be 
calculated or modeled. These correlations correct the 
laminar chemical reaction rate for turbulence. All of these 
correlations involve modeling assumptions . The pdf method 
does not include these assumptions. This is why pdf 
modeling is promoted as exactly modeling chemical reaction. 
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The terms on the right side of the pdf evolution 
equation must be modeled. The turbulent convection term is 
modeled by gradient diffusion: 

-<»?!**>£ - D c aj, (4 ' 4) 

where D t is the turbulent diffusion coefficient, which is 
related to the eddy viscosity by the turbulent Schmidt 
number. The next term that needs modeling is the molecular 
mixing term. A proper molecular mixing model causes the 
scalar variance to decrease with time while the mean remains 
constant. Some molecular mixing models are: the 
coalescence/dispersion model attributed to Curl 71 , modified 
Curl model of Janicka et al. 72 , and the quasi-Gaussian 
model of Dopazo 73 . The general form of the 

coalescence/dispersion model is: 

pi) £ < <«*> !♦*>*> “-T // i3r( * ,) * ( * ,/)r( *'^ t'W <**'<**" - p 

i=l J=1 1 

■ M(P) 

( 4 . 5 ) 

where T is the probability that particles will undergo a 
transition in concentration between one particle with 
concentration of Y and another particle with concentration 
\|/» to produce two particles with concentrations and \j/'+\(/ 

\|/ respectively. In Janicka et al.'s paper, a general 
equation for transition probability is developed using the 
concentration mixing equations developed by Curl for droplet 
mixing. Using the simplest approximation, the transition 
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probability is: 


r(V f <I» // N' i ) =i 


(ti - Vi) 

o 


fortyUVi^'l or 


(4.6) 


otherwise, 

This approximation assigns equal transition probabilities 
for each point in the range \|f' to \|T , The time scale x for 
this process is taken as the integral time scale, T=k/e. 

Traditional finite difference techniques are not 
feasible for solving this pdf evolution equation as the pdf 
is a function of a large number of variables (five species 
plus temperature or enthalpy) . Instead, a Monte Carlo 
method will be used to solve for the pdf evolution. 

4.1 - The Monte Carlo Solution Method 

The Monte Carlo method is easy to implement, even if 
the number of variables is large. The evolution of the pdf 
entails the movement of particles in physical space as well 
as the scalar space (\j/ space) . The movement of the 
particles is governed by the pdf evolution equation. The 
solution to the pdf evolution equation is a finite set of 
particles or events for each grid cell. The actual pdf is a 
continuous function. Increasing the number of particles 
does not result in a continuous solution. Another problem 
with Monte Carlo is the level of solution accuracy. The 
ensemble average a solution particle property (species 
concentration, temperature, enthalpy or density) for a grid 
cell is an estimate of the true mean since the expected 
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value of the ensemble average is the mean. It can be shown 
that the solution accuracy is proportional to the square 
root of the number of particles . Thus to get an order of 
magnitude increase in accuracy the number of particles must 
be increased by two orders of magnitude. 

Mathematically, the particles for the Monte Carlo pdf 
are represented by a set of delta functions . The 
distribution function is the integral of the probability 
density function. The Monte Carlo distribution function is 
a sum of heavy-side functions. An example from Pope's 1985 
paper 74 is shown in figure 4.1. A Gaussian pdf and a set 
of 20 normally-distributed random numbers are compared. The 
mean is 0.5 and the standard deviation is 0.01 for both 
probability density functions. It can be seen that the 
Monte Carlo distribution function approximates the 
continuous distribution function in a step-wise manner. The 
delta functions are represented by finite spikes. The 
actual height of the spikes is infinite and the integral of 
each spike is 1/20. 

Using fractional time-steps, the pdf solution is 
divided into 3 phases. The time derivative is first 
discretized as: 

(jP ° +1 —p n ) / At (4 - 7) 

then the pdf evolution equation can be written as: 
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Figure 4.1 Comparison of continuous and discrete pdf's and 
distribution functions 
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P n+1 = [1 - AtpO s a a - A tV d 0i (* x , . . . + A td a D t d a + A tM] P D 

( 4 . 8 ) 

Using approximate factorization, the above equation is 
recast as 

P“ + Ml-At£ . . . ,*,) > (1+AtM) {l-AtpVd a +±td a D t d a ) P n 

i=i 


«(l+AfcO (1+AtM) (1+AtP) P D ' ' 

where C denotes the convection/dif fusion operator, M is the 
molecular mixing, and R is the chemical reactions. Using 
the above expression the different processes may be 
performed consecutively: 

(1) convection/diffusion between grid cells: 


P* = (1+AtO P D , 

(2) molecular mixing within each grid cell: 


( 4 . 10 ) 


P** = (1+AtM) P*, 


( 4 . 11 ) 


(3) chemical reaction in each grid cell: 


P " +1 = (l+At P)P**. 

In the Monte Carlo pdf calculation, each particle has 
its own species mass fraction and temperature -enthalpy. 

Each particles has the same mass. Particle position within 
gridcells is not modeled. In the Monte Carlo simulation a 
continuous pdf is replaced by N x m delta functions, 
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!>*(*„*, *„> - -sE 8 (t!-*!”’ ) *« ) ~» ) (4.13) 

"j>1 

each product of the m delta functions representing one event 
of an ensemble of N sample events. The pdf is a function of 
position and time. 

The Monte Carlo pdf formulation/solution is initialized 
with a SIMPLE method solution. The pdf solution technique 
could be started with an unconverged flowfield, but, the 
time needed to get a converged solution would be much 
higher. Each of the initial particles in a grid cell have 
identical properties. Inlet particle properties are set to 
the inlet boundary conditions. The first step in the 
iteration technique for the pdf solution is calculating the 
exchange of particles between grid cells by convection and 
diffusion . 

4.1.1 - Convection/Diffusion phase 

Replacing the continuous pdf with the Monte Carlo pdf, 
the convection/diffusion step of the evolution equation can 
be written as: 


p3 t P* + poa.p* = dj (D t d a P*) (4.14) 

re-writing this for the x-y coordinate system: 


app* 

at 




jL_a_ 

y Sy 


yp VP* - yT p 


dP* 

dy J 


= 0 


(4.15) 


The convection/diffusion step does not have any source 
terms. The SIMPLE scheme solves the transport of species 
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semi-implicitly . That is, very large time steps can be 
taken. Solution information travels very quickly in this 
system. The Monte Carlo pdf scheme uses the evolution of 
particles to calculate scalar properties such as 
concentration, temperature and density by ensemble 
averaging. The Monte Carlo scheme does not implicitly 
calculate particles or particle properties. Particle 
properties are calculated from existing neighbor particles 
in time and space. Thus, it is incorrect to use implicit 
differencing in Monte Carlo simulations. Explicit 
differencing must be used, which limits the time step that 
may be taken. It should be obvious that the time step must 
small enough that any particle does not travel further than 
one grid cell in any direction. In the SIMPLE scheme the 
weighing of neighboring grid cell properties is based on 
coefficients calculated from the convection and diffusion 
using an implicit formulation. The coefficients used in the 
Monte Carlo solver must be figured using explicit formulas 
as the exchange of particles is explicit. For example, 
assuming 1-D diffusive heat transfer we have as a general 
formula : 

a P T p =a s [fT B + (1 -f) Tj] +a w [fT w + (1 -f) r£] + [a$- (1-f) a B - (1-f) a v ] T P 

( 4 . 16 ) 


where 
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a p = fa s + fa w + a% (4.17) 

and a E and a w are calculated as before. The letter f is the 
weight of implicit differencing. Generalizing this for 2- 
dimensional convection and diffusion with f=0 for explicit 
differencing: 

(4.18) 

where 


„o _ pc ax Ay 
3 - ■ *» - At 


(4.19) 


thus, the time step must be small enough so that the 
weighing for the previous time step is non-negative (which 
introduces non-physical solutions) . A maximum step is 
calculated for each iteration by finding the cell with the 
smallest allowable time step and using that value. 

Dividing the above equation by a p and using the Monte 
Carlo pdf with i,j notation, the convection/diffusion step 
is : 


Pl.j = aP i*i,j+i + PPi-i.j-i + 1 + ep i.j 




In this step the new Monte Carlo pdf is made of a fixed 
number of particles. This is a combination of particles 
from neighboring cells and the particles from the same cell 
(particles from the previous time step) . The particles from 
each cell are randomly chosen. Only whole numbers of 
particles are exchanged. In the conversion from a floating 
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point number into an integer number some fraction is lost. 

The coding is such that the remainder is added in at random 
iterations so that the exchange process is more accurate. 

If the remainder were not randomly added in, it is thought 
that the iteration process may become cyclic, at least with 
a small number of particles per cell. The time step is 
reduced so that the possibility of extra convection or 
diffusion exchange particles does not cause the number of 
particles at the previous time to go non-negative (ie. 
exchange more particles than you actually have) . Of course, 
using high numbers of particles per cell cuts down on the 
truncation error. The Monte Carlo scheme does not keep 
track of where the particles are within each cell . The 
solution quantities to a SIMPLE method calculation (as most 
CFD solutions) represent the means and sometimes the modeled 
variance for individual cells. The variance represents the 
degree of mixing within a grid cell . 

4.1.2 - Molecular Mixing Models 

In the second fractional time step, molecular mixing is 
performed using the continuous mixing model of Hsu. The 
modified Curl model has non-gaussian higher moments. Curl s 
model assumes complete mixing between chosen mixing pairs. 

In the modified method of Curl, the Monte Carlo simulation 
randomly chooses pairs of particles based on the 
following formula: 
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N— = 0.5 = 0.5-^-N, 

mx Cx Ck 


(4.21) 


where C is 6.0. The pairs are allowed to mix as: 

* n (t+ftt) = Ay a (t) + (1-A)* a (t) (4.22a) 

y m (t+6t) * AV n (t) + (l-A)$ Ia (t) (4.22b) 

where A = 0.5£, with £ a random variable uniformly- 
distributed on the interval [0,1]. The remaining N - 2N mx 
particles remain unchanged: 

i|r B (t+$t) * f n (t) (4.22c) 

A plot of the mixing process in composition space is shown 
in figure 4.2. The value of A used in this plot is 1/6. 
This model is discontinuous as the sample particles change 
properties abruptly regardless of how small the time 
interval is made. Reducing the time interval reduces the 
number of particle pairs chosen for mixing, but, does not 
affect the total property change for the chosen particles. 
A. Hsu and J.Y. Chen 75 modified this process to become 
continuous. In this model all particles are chosen to 
undergo mixing and the time interval is incorporated into 
the mixing level term A. Thus, the jump in physical 
properties mixing pair particles undergo decreases with the 
time step, and thus, the process becomes continuous. There 
are N/2 particle mixing pairs for each cell. A is now 
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Figure 4.2 Plot of Monte Carlo Mixing process 


defined as: 

A = C>i±j±, ( 4 - 23 ) 

where C 7 = 2.0. Hsu showed that this method resulted in 
more appropriate higher order moments, specifically flatness 
and superskewness in his calculations. In a 1990 paper, 

Chen and Kollman 76 , use the modified Curl method for a two- 
step hydrogen combustor calculation and conclude that the 
model did predict cross correlation properties reasonably 
well for turbulent shear flows, despite the problems with 
non-Gaussian character in the higher moments. 
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Other mixing models have been suggested and employed. 
Pope 77 suggested that the sampling probability of particles 
should be biased in regards to the elapsed time since the 
last mixing event. This entails an additional scalar for 
the elapsed time since the last mixing event for each 
particle and the use of an appropriate weighing function for 
age. Another mixing model is the relax-to-mean model. In 
this model, the cell mean is calculated, then each 
individual particle's variance from the mean is reduced 
exponentially depending on turbulent mixing and time step. 
The formula for this is: 

i|r n ( t+6 t) =i|f+ (i|r n ( t) -i[r) *exp (-0 . SC'tS t) (4.24) 

This mixing model was used in A. Norris and A. Hsu's 78 
Monte Carlo pdf calculation of Masri et al 79 CH 4 /CO/H 2 /N 2 - 
air pilot stabilized jet flame. A somewhat similar model 
which uses a linear relaxation is the IEM model developed by 
Villermaux , 80 and by Yamazaki and Ichigawa 81 . The formula 
for this is: 


dYi 

at 


+ u 


dYj 
3 dxj 



+ o i 


(4.25) 


The model name comes from exchanges by Interaction with the 
Mean Value. The term t iem is an exchange time. This model 
treats combustion simultaneously with mixing. This model 
has been used by Correa and Braaten 82 to test reduced 
chemistry mechanisms against full chemistry mechanisms for 
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methane-air combustion in a Partially Stirred Reactor (PaSR) 
model. This model was chosen in preference to pair exchange 
models because it is easier to employ on parallel computers. 
The simulations take many hours to perform on a single 
processor. Doing calculations in parallel, that is, 
dividing the numerical work between multiple machines will 
result in greatly improved turn-around time. The PaSR 
represents a step into pdf modeling which can use full 
laminar chemistry for a single grid cell or reactor. 

The PaSR pdf calculation admits a fixed number of inlet 
particles per mixing step. At the same time, the same 
number of particles is taken from the PaSR. The PaSR has 
two other simpler reactor models as mixing limits. With no 
mixing the reactor is simple plug flow. The other limit is 
infinite mixing which corresponds to a fully stirred 
reactor. Actual multi-dimensional flowfield calculations 
can be thought of as a collections of Partially stirred 
reactors. The exchange of particles between grid cells 
amounts to exchange of particles between neighboring 
reactors. The mixing frequency inside each reactor is 
governed by local turbulence mixing. The results of Correa 
and Braaten' s PaSR calculations varied with the frequency of 
mixing. Specifying low mixing frequencies resulted in 
blowoff for a reactor with separate fuel and oxidant inlets. 
Even with steady-state combustion, concentrations, 
particularly NO and CO, varied with mixing frequency even 



using full-up chemistry. As reducing CO and NO emissions is 
of prime concern in engine design, this is somewhat 
troubling. Taking individual gas samples amounts to 
measuring an average. This also disturbs the flow, which 
affects accuracy. Laser Raman Spectroscopy can take nearly 
instantaneous concentration measurements without disturbing 
the flow, but, the measurements must undergo substantial 
corrections. Thus, it may be some time before data are 
calibrated well enough to fully validate the pdf method. 

In a later paper, Sanjay Correa 83 compared the IEM 
mixing model with Curl's model and the modified Curl's model 
for premixed CO/H 2 combustion in a PaSR. At a frequency of 
316 Hz, the results were similar. From figures 4-6 of his 
paper, the difference in temperature between IEM and the 
pair exchange models was on the order of 50 degrees. The 
difference in predicted CO and OH concentrations was close 
to a factor of two. Oxygen predictions were not affected by 
the mixing model at this frequency. At a frequency of 1000 
Hz, the pair exchange models and the IEM model produced 
nearly identical results . As Correa and Braaten predicted 
in-flame mixing frequencies of 1000 Hz and larger, Correa 
stated that the choice of mixing models was not critical in 
the distributed reaction regime of lean premixed combustion, 
as long as the turbulent mixing frequencies were above 1000 
Hz. This thesis will employ benchmark data from research 
combustors. These combustors operate at lower temperatures 
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and less intense combustion than practical combustors. 

Thus, the I EM model is probably not appropriate to use for 
this work. Research combustor data is used in this work as 
experimental data from practical combustors is sparse. 

Combustion was also effected in a Coalescence- 
Dispersion Model used for infinite reaction cited by J.Y. 

Chen and W. Kollman 84 . If chemistry is very fast compared 
to mixing, as in the flamelet combustion modeling, 
combustion occurs within a very limited mixture range close 
to stoichiometric. If combustion isn't treated 
simultaneously with molecular mixing, combustion can be 
vastly under predicted. In scalar space, if only the final 
mixing concentration is examined, particles can pass through 
a combustible mixture during the mixing phase without 
reacting. This problem does not occur with distributed 
combustion which is considered for this thesis. 

4.1.3 - Combustion phase 

Combustion is proceeds after the exchange and mixing of 
particles. The pdf method does not have to employ 
corrections to reaction rates due to turbulence. The 
unmodified laminar reaction rate is used in the pdf 
evolution term. Each particle has its own species 
concentration and enthalpy. Each particle then undergoes 
chemical reaction without having to worry about modeled 
turbulence correlations, turbulent flame speed or other ad- 
hoc turbulent combustion rate. The laminar reaction rate is 
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the actual reaction rate. However, as mentioned previously, 
there is much controversy about true reaction rates. Even 
laminar combustion involves tens of reactions with diverse 
time constants. Including "exact" kinetics into a pdf 
scheme would result in 100,000 "exact" solutions per 
iteration for a rough grid of 50 by 80 with 250 particles 
per grid cell. Semi— implicit turbulent combustion 
calculations take thousands of iterations to converge a few 
orders of magnitude. Explicit pdf calculations require at 
least ten thousand iterations for elliptic flows. The 
current use of pdf methods must inherently include 
simplification combustion modeling to give realistic 
calculation times. Some of these approximations have been 
already mentioned. If the number of scalars is kept low it 
is possible to construct thermodynamic state tables. The 
table is then interpolated to find thermodynamic data during 
the iteration process. For instance, J. Y. Chen et al. 85 
have constructed state tables for methane combustion using 4 
and 5 scalar reaction mechanisms and two time increments. 
Reaction submodels incorporated in the Lewis pdf model have 
included infinite reaction rate, a two step global reaction 
model for hydrogen by Chinitz, an Intrinsic Low-Dimensional 
Manifold (ILDM) method of Maas and Pope 86 . The ILDM method 
uses dynamical systems approach to reduce the number of 
degrees of freedom to simplify the chemical kinetics system 
and construct a state look-up table. The method was 


HIT!" 



83 


successful in the calculation of Norris and Hsu [78] . The 
ILDM method for heavier, more complex fuels more difficult 
and is currently being worked on. Also, the table look-up 
method for higher molecular weight, more complex fuels 
requires more reaction steps. This is beyond current 
computer memory limits. A single finite-rate reaction 
submodel is used in this work. 

In 1990, Chen and Kollman 87 compared simplified 
chemical reaction models for hydrogen-air combustion. At 
conditions near equilibrium, the chemistry models all 
produced similar results. Only one of the models was 
reported to be able to predict flame blowout, but, the 
blowout velocity was incorrect. They concluded that further 
improvements in the turbulent mixing model and the numerical 
scheme were required to obtain quantitatively correct 
predictions. Also needed are data from simultaneous 
measurements of radical species and mixture fraction. 

4.2 - Further Development of Composite PDF Method 

One contribution of this work is further development of 
the composition pdf and testing the model for finite-rate 
combustion in recirculating flows. The composition pdf has 
been primarily used for parabolic or one-way flows with 
modifications of equilibrium type combustion models. 
Equilibrium or fast kinetics combustion models should be 
used where the rate of mixing is much smaller than the 

These models are valid for low molecular 


combustion rate. 
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weight fuels as hydrogen and methane at high temperatures. 
Practical combustors use higher molecular weight fuels. 
Burning is distributed in practical combustors. There is 
significant overlap in reacting species, particularly at low 
temperature . 

The SIMPLE method of Patankar and its variants is the 
most widely used solution method for practical combustor 
flows. The composition pdf method is combined with the 
SIMPLE method. Momentum, overall mass conservation, 
turbulent kinetic energy and dissipation are solved using 
the SIMPLE method. The composition pdf model used in this 
work was initially developed by Hsu and co-workers at NASA 
Lewis Research Center. The first release of this 
composition pdf model was with the Lewis RPLUS 2-D combustor 
code 88 . RPLUS is used to predict high speed combustor 
flows. As Mach number is lowered the solution convergence 
rate deteriorates. It is an extension of the Lower-Upper 
Symmetric Solver used by Yoon and Jameson 89 to predict 
species composition. The first dissemination of the 
development model was at a workshop in 1993. The model 
released at the workshop was set up to solve the flow 
configuration of Burrows and Kurkov' s supersonic hydrogen 
combustor 90 . The model was set up to use a global five 
species, two step reaction submodel originally proposed by 
Rogers and Chinitz 91 for hydrogen-oxygen combustion. 

Many changes were made to the original composite pdf 
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model. The implementation of the finite rate combustion 
model didn't calculate steady state combustion using 
suggested reaction numbers for methane. A different species 
composition solver was implemented. The species 
concentration solver uses the Newton-Raphson technique. In 
the original coding, an incorrect implementation of upwind 
differencing was used to calculate the number of particles 
to exchange between cells. Pure upwinding is regarded as 
highly diffusive in CFD . So as to fairly compare the 
predictive capability of the combination pdf model against 
the baseline combustion model, hybrid differencing was 
implemented in the composition pdf for the calculation of 
exchange coefficients. The thermodynamics of the base pdf 
model were highly simplified. Constant C p was assumed. C p 
changes markedly over the range of combustor temperatures. 
Predictions of exit temperature were adjusted in development 
model by adjusting C p . The baseline combustion model used 
in this work has variable C p . Variable C p was also 
incorporated into the thermodynamics of the composition pdf 
model. This necessitated the use of enthalpy rather than 
temperature as one of the pdf variables. The development 
combination pdf model used a combination of temperature 
solver in Monte Carlo solver and enthalpy solver in RPLUS 
code to solve for temperature. The model implemented here 
properly solves for temperature completely in the pdf 


solver . 
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Trial calculations are done using the SIMPLE method as 
it is more economical than the hybrid pdf method. The best 
solution is used to begin hybrid pdf solution and is also 
used to compare predictions and measurements. 
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Chapter 5 

COMBUSTOR MODEL VALIDATION 
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Experimental data used to evaluate modeling should be 
selected based on the quality and quantity of experimental 
data, especially initial or boundary condition data. 
Experimental data should include measurement of multiple 
chemical species, velocity, temperature and turbulence 
quantities. If comparison is performed with limited data, 
the quality of a CFD solution may not be adequately 
evaluated. CFD solutions are usually repeated changing 
initial conditions, gridding, turbulence modeling, etc. 
until an optimum comparison is given. With limited data, an 
optimum solution may be found for the wrong reason. A good 
prediction of one species does not necessarily mean that all 
other species are equally well predicted. Similarly, a good 
prediction of one velocity profile does not mean that other 
velocities profiles are equally as good. Good predictions 
at one station in a combustor does not mean that all 
stations are of the same quality. 

Suitable and accurate data must be used to ascertain 
the appropriateness and accuracy of numerical predictions. 
Past measurement techniques have used hot wires and pressure 
probes to determine velocities. Over the past decade, 
optical measurements of velocity have become standard. 
Measurements require some kind of access. Probes can be 
handled by access from the rear of the combustor. Lasers 
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demand optical access through the combustor. This is 
usually handled by installing heat resistant windows in the 
combustor. The windows degrade faster with increasing 
temperature, heat transfer and pressure. Thus, almost all 
benchmark quality data for combustors have assembled using 
research-type combustors. These combustors incorporate the 
major features of practical commercial gas turbine 
combustors, but, their designs are vastly simplified. 

Running an actual combustor is very expensive and very 
difficult to instrument. Most research type combustors are 
axisymmetric, except for the possible addition of dilution 
jets. The dilution jets are required to cool the combustion 
liner and lower the overall equivalence ratio of the 
combustor. The solution of combustors with dilution jets 
requires 3-D models. Combustors flows undergo substantial 
change in density. The flow through a combustor without the 
flow expansion caused by combustion can be substantially 
different than with combustion. Measurements of some cold 
combustor flows have shown reverse flow at the exit of the 
combustors. Most CFD models assume one-way flow at the 
solution exit plane. One way exit flow may be enforced by 
installing flow restrictions at the combustor exit. 

5.1 - Combustor Model Validation Data 

In 1985 an evaluation of data for parabolic turbulent 
reacting flows was published. The editors were Strahle and 
Lekoudis 92 . The work was sponsored by the Air Force Office 
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of Scientific Research. The purpose was to conduct a 
program similar to that run at Stanford in 1968 93 , except 
on turbulent reacting flow as opposed to turbulent boundary 
layers. A committee was set up to help in the evaluation of 
the data. At the final committee meeting it was recommended 
that a computational effort not be initiated at the time. 

One reason was that most combustion theories are application 
specific and cannot be readily used for flows of different 
character or chemistry from those for which they were 
developed. Some problems noted in the data were: lack of 

completeness, incomplete specification of initial and 
boundary conditions and excessive experimental uncertainty. 
This work did not include data for complex reacting 

turbulent recirculating flows. 

One of the major problems with early combustion data is 
the limited number of experimental quantities measured. 

Many of the experiments were performed to characterize only 
parametric variations of a combustor type, not to validate 
CFD models. In a paper by Baker et al. 94 , only velocities 
and stresses are reported. In another paper by Lockwood et 
al. 95 , only flame length and mixture fraction were 
reported. McDannel, Peterson and Samuelson 96 gave species 
and temperature measurements for a premixed flow. The flame 
was stabilized by a recirculation zone formed by a high 
velocity jet opposed to the direction of bulk flow. 
Velocities were varied in a the study. El Banhawy, 
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Sivasegaram and Whitelaw 97 took measurements of 
temperature, velocity and species concentration for a 
premixed combustor in 1983. The measurements are presented 
in contour form. In premixed combustors the fuel and 
oxidant are combined upstream of the combustor inlet. In 
this configuration, combustion is anchored by a hot 
recirculation zone caused by a sudden area increase. This 
geometry is known as a step combustor. Two different step 
sizes were studied. Premixed combustors aren't used in 
practical gas turbine engines, due to stability problems. 

The concept has been under study for quite some time because 
the concept reduces engine emissions. Smith, Giel and 
Catalano 98 measured hydrogen-air combustion in an 
axisymmetric recirculating combustor flow. A lean mixture 
was used. Due to its low molecular weight, hydrogen 
diffuses much more rapidly than oxygen or water. 

Preferential diffusivity wasn't included in the numerical 
model, so this flow was not calculated. There is a wide 
disparity in the inlet velocities, which causes the 
recirculation zone. The hydrogen has a bulk velocity of 
.914 m/s, while the air inlet velocity was 102 m/s. 
Temperature and species concentration for the reacting case 
was measured by laser raman spectroscopy. However, not 
enough measurements were taken at each location for high 
confidence levels. The temperature profiles are quite rough 
for the first three measurement stations. 
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Lightman et al." performed measurements on a bluff 
body combustor at Wright-Patterson Aero Propulsion 
Laboratory. Measurements were for combusting and non 
combusting flows. They reported that most of the features 
of the cold flowfield also were present in the combusting 
flow. However, in a later paper 100 , high speed pictures of 
the combustor flow revealed the existence of large scale 
structures. By examining laser Doppler anemometer 
measurements and measurements of luminosity it was concluded 
that the flow is comprised of non-luminous regions and flame 
turbules. The combusting flow was unsteady. Schefer et 
al . 101 also saw coherent structures in their study of 
another combustor flow stabilized by a bluff body. The flow 
was characterized by large scale entrainment of air into the 
recirculation zone. Pdf's of Velocity and velocity 
correlation were measured. Bimodal pdf's were observed in 
regions of high shear in the recirculating zone boundaries 
and in the downstream stagnation region. This bimodality 
was associated with the alternate passage of large-scale 
turbulent structures consisting of unmixed fuel and air 
through the measurement volume. In 1987 , Sivasegaram and 
Whitelaw 102 published a paper studying oscillations in 
axisymmetric dump combustors. Another study was published 
in 1989 by Sivasegaram et al . 103 Some unsteady combustion 
predictions have been done, but, the practicality of the 
work is somewhat controversial. 
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Samuelson and his students have produced a number of 
papers regarding the design and experimental measurements of 
axisymmetric combustors. The most complete 2-D database is 
an University of California , Irvine, Combustion Laboratory 
paper by R. Charles 104 . A parametric series of 
measurements were made by modifying the combustor inlet 
conditions including the fuel nozzle. Unfortunately, only 
axial and tangential velocity measurements are reported. 

The radial velocity is unknown, although it may be inferred 
from continuity for non-combusting flows. In general, the 
number measurements taken in the radial direction number 
around ten. This number is too small to accurately validate 
or benchmark combustor models. Axial velocities measured 1 
cm from the inlet show considerable backflow. Azimuthal 
velocity generally shows a peak near the outside wall at 
this axial location. Peak azimuthal velocity at other axial 
locations is much closer to the combustor centerline. Inlet 
data from this configuration would have to be set up based 
on flow geometry, bulk flow and experience of the combustion 
modeler. 

J.C. Pan 105 did a parametric study of turbulent 
confined premixed flames. The flame was stabilized by 
conical bluff bodies. In the study, the size and angle of 
the conical bluff bodies, the equivalence ratio, the inlet 
velocity and turbulence were varied. The data given in the 
reference were not in tabular form. Velocity and 
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temperature pdf's were given, velocity vector diagrams, 
contours of velocity and reynolds stresses and centerline 
profiles of velocity and temperature. Some cases used very 
lean mixtures close to the blow-out limit. The data may be 
good for validating flamelet combustion models. 

Some three dimensional combustor flows have also been 
measured. Heitor and Whitelaw 106 measured velocity, 
temperature and species in a combustor having dilution jet 
holes. The combustor exit is constricted to one side. The 
fuel was gaseous propane. The scalar field was dominated by 
the presence of the dilution jets. Air to fuel ratio was 
varied. Combustion efficiency was low for some of the runs. 
In a 1990 paper, Bicen, Tse, and Whitelaw 107 reported up to 
98 % combustion efficiency for a similar can-type combustor. 
The primary zone of the combustor was fuel rich. 
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5.2 - Spray Model Validation 

The baseline combustion model was tested in various 
stages to benchmark its predictive capability. The first 
stage was predicting an isothermal open combustor flow. The 
spray vaporization model was then tested. Finally, 
calculations of spray combustion were performed. 

5.2.1 Isothermal Swirling Flow 

The isothermal swirling combustor flow of D. Bulzan 108 
was used to validate the baseline combustor model. This 
combustor is part of an on-going effort to provide experi- 
mental data to validate spray combustion models. The exper- 
imental configuration is a Parker Hanifin research simplex 
air-assist atomizer surrounded by a co-flowing air stream. 
The configuration used here is unconfined, that is, it does 
not have an outer wall. This was done to simplify taking 
experimental measurements . A schematic drawing of the 
combustor is shown in fig. 5.1. Most unconfined or free- jet 
combustors are notorious for poor mixing, and require sig- 
nificant axial distances for complete mixing. The atomizer 
and the co-flow air had 45 degree swirl in the same direc- 
tion. This swirl causes the formation of a large recircula- 
tion zone which promotes rapid mixing and stabilizes the 
flame for the reacting case. The recirculation bubble is 
longer than the laser measurement system can traverse. The 
model was able to handle this negative outflow condition by 
suitable treatment of the exit boundary conditions as sug- 
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gested by P. George 
Huang. Normally, exit 
boundary conditions 
assume zero-gradient or 
one-way flow where the 
flow is not signifi- 
cantly changing with 
axial distance. The 
exit velocity is un- 
known, but, is needed 
in the calculation of 
pressure. The solution 
is part of the itera- 
tive process. The exit 
velocity is set to the 
last calculated ups- 


— Atomizing 
air swirler 



-Coflow 


Figure 5.1 Schematic drawing of the 
combustor. Dimensions in mm. 


tream velocity plus a 

velocity correction due to mass conservation. Thus pressure 
for the exit computation cell can be calculated. In the 
process of calculating the upstream velocity, the downstream 
contribution to the exit computational cell is set to zero 
without adding any significant error, since the downstream 
contribution is either zero or very small. For the case 


where the exit velocity is negative for a computational 
cell, the exit neighbor coefficient is the upstream coeffi- 
cient and the contribution is very significant and cannot be 
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ignored. This is actually an unspecified inlet condition, 
which must be calculated. This is handled in the iteration 
process by setting the boundary cell velocity equal to the 
velocity of the downstream adjacent cell, which is calculat- 
ed. Then the upstream or boundary neighbor coefficient is 
added to the non-constant part of the source term. Then the 
neighbor coefficient is set to zero (the boundary coeffi- 
cient isn't included in the sum of neighbor coefficients). 

To summarize, first neighbor coefficients are calculated in 
the normal way. Then the (negative) exit velocity coeffi- 
cient is added to the non-constant part of the source term. 
Then the exit velocity is treated as if it were zero in the 
difference equation. Pressure is still calculated the same 
as before. The mass correction at the exit boundary is only 
made where the exit velocity is positive. 

The flowfield in the inlet region is quite complicated. 
The combination of the blunt body atomizer with the highly 
swirling co-annular flow produces two recirculation zones. 
The first has already been mentioned. The second is a small 
toroidal vortex adjacent to the blunt atomizer between the 
outward swirling flow and the reverse flow along the center- 
line. The unconfined combustor entrains surrounding ambient 
air. A velocity-vector diagram from a calculation is shown 
in figure 5.2. The centerline is on the right side of the 
figure. The flow is from bottom to top. The small vortex 
near the bottom center-line is somewhat larger than what was 
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experimentally measured. Calculation of the flow-field is 
difficult. One problem is reverse flow at the first mea- 
surement location. The current model can handle moderate 
amounts of reverse flow at the combustor inlet, as it uses a 
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pressure based algorithm. Another problem is that the ini- 
tial flow is rapidly expanding into an open environment. 
Entraining flow is difficult to accurately predict using 
pressure based algorithms. The strong swirling flow, which 
produces the very large recirculation bubble, could not be 
contained by boundary conditions intended for entraining 
free- jet flow. Using wall-type boundary conditions did 
contain the flow, but, the size of the calculated recircula- 
tion zone depended on placement of the non-physical wall. 
Using a larger radial calculation domain produced a larger 
recirculation bubble, with most of the positive axial flow 
occurring close to the wall . The experimental measurements 
show an axial velocity peak which gradually moves outward. 

A solution to this was to specify velocities on the outer 
radial boundary. Dr. Bulzan measured some velocities at 
radial distances of 143 mm and 168 mm from the combustor 
centerline. This data was obtained directly from Dr. Bul- 
zan. The measurements at 143 mm were used for this calcula- 
tion. At this radius, flow was entrained to an axial dis- 
tance of about 130 mm. Measured values of velocity at 5 mm 
downstream of the atomizer were used as initial conditions 
for the calculation. Predicted flow fields varied greatly 
by using different initial values of turbulent kinetic 
energy and dissipation. Normally, the k-£ turbulence model 
is believed to produce only small change due to rather large 
change in initial turbulence quantities. Typically, the use 
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of poor initial turbulence quantities is overcome within a 
few gridpoints. Thus, the k-t model is perceived as a good 
model to use as predictions are not unusually dependent on 
initial dissipation or length scale, which are usually 
unknown. Unfortunately, the initial structure of the flow- 
field is quite complicated, is unconfined, and involves very 
high shear. The atomizer causes large velocity gradients. 
This causes a huge generation of kinetic energy and rapid 
dissipation. Using the full measured turbulence energy and 
a length scale of 30 % of the diameter of the atomizer and 
co- flowing air passage smeared out the small secondary 
vortex and blocked the reverse flow along the centerline. 
Also, the predicted central recirculation zone was over a 
factor of two wider than measured at the calculation exit. 

Another calculation was done using peak turbulent 
kinetic energy around an order of magnitude smaller than 
measured and specifying a length scale of 30 % of the width 
of the total calculation domain. This resulted in very good 
predictions of the velocity field. Evidently, the low ini- 
tial values allowed lower initial dissipation and which 
later built up enough to predict the correct amount of 
spreading . 

Due to convergence problems, initial swirl was slowly 
ramped up to the full value. The total number of iterations 
was 6000. Calculated versus measured axial velocities are 
shown in figure 5.3. The comparisons are done at axial 
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distances of 20, 50, and 100 mm. The axial velocity predic- 
tions show quite good agreement with the measurements. The 
velocity peaks are slightly underpredicted. Figure 5.4 
shows comparisons of radial velocity. Radial velocities are 
underpredicted . 

Figure 5.5 shows comparisons of tangential velocity. 
Similar to axial velocity profiles, the calculated tangen- 
tial velocity profiles are very good. Peak velocities are 
only slightly underpredicted. Plots of turbulent kinetic 
energy are shown in figure 5.6. The comparison is poor for 
the lower half of the profile at 20 mm, where the turbulent 
kinetic energy is vastly underpredicted. This part of the 
flow represents shear between the second vortex and the 
backflow along the centerline from the large recirculation 
zone. The local axial velocity peak is also underpredicted 
in this shear region. 

Predicted profiles at other locations are very good. 
These calculations show that juggling initial turbulence 
values has a significant effect on computational results for 
this flow. As reducing initial turbulence values reduced 
mixing and produced good computational results, a modifica- 
tion to the turbulence model was considered. The k-e model 
assumes isotropic turbulence and is not regarded as optimal 
where swirl, body forces, or curvature exist. In the past, 
some researchers have resorted to a curvature corrections to 
improve predictions for strongly swirling flows 109 . 
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Figure 5.4 Comparisons of radial velocity 
at 20, 50, and 100 mm. 
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Several modifications were tested by Sloan et . al. no for 
thirteen cases. Unfortunately, the lack of initial condi- 
tions prevented drawing absolute conclusions about the 
improvements . This work does not employ curvature or Rich- 
ardson number corrections. Instead, a simpler method of 
changing the turbulent viscosity coefficient based on shear 
was tried. Recently, a modification for high shear was 
proposed and tested by Shih et. al. 111 . This development 
is based on invariant theory in continuum mechanics. The 
modification is done so that negative normal stresses are 
not calculated in any situation of rapid distortion. This 
involves damping of the coefficient for turbulent viscosity, 
which is normally treated as a constant. The equations used 
are: 


v =c-£i 
' f e 

(5.1) 

C = 2/3 

(5.2) 

“ Ai+Tl 


■q = Sk/e 

(5.3) 


(5.4) 


(5.5) 


The predicted viscosity coefficient ranges from a high of 
0.121 to a low of 0.0178 at the edge of the co-flowing air. 
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The results of this calculation are shown in figures 5.7-10. 
The axial velocity profile at 20 mm is similar to the 
results with the regular k-e model, but, the 50 and 100 mm 
profiles do not show as much radial spreading. The radial 
velocity profiles are slightly reduced which accounts for 
the lower spreading rate. This is also reflected in the 
tangential velocity profiles. The peak turbulent kinetic 
energy for the 20 mm axial location is improved, but, at all 
other locations predictions are worse with much less radial 
spreading. It appears that modification of the turbulence 
can be made to improve predictions in one area, but, unfor- 
tunately reduce accuracy in other areas. 

These calculations were done with hybrid differencing 
which is known to be overly diffusive, especially for calcu- 
lation of round free jets. Higher order differencing should 
be used before judgements about the validity of various tur- 
bulence modeling improvements. 
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Figure 5.7 Axial velocity profiles at 20, 
50, and 100 mm with modified turbulence 
model 
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Figure 5.8 Radial velocity profiles at 20 
50, and 100 mm. 
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5.2.2 Evaporating Spray Model Validation 

The methanol vaporization data of McDonell et. al. 112 
was used to validate the spray model. McDonell' s data were 
specifically taken to validate spray models. The experimen- 
tal apparatus consists of a research atomizer injecting 
methanol spray downward into a 457 mm square duct with co- 
flowing air stream at an approximate bulk velocity of 1.0 
m/s. Measurements were made of axial and radial gas-phase 
and droplet velocities, methanol vapor concentration, and 
droplet fluxes. No reverse flow was measured in this non- 
combusting case. The droplet data were measured by Phase 
Doppler Interferometry (PDI) . The spray droplet data were 
partitioned into 10 different size groups due to limitations 
of the instrument. Data at an axial distance of 7.5 mm were 
used as initial conditions for the calculation. Gas phase 
velocities were measured to a radial distance of 88.0 mm. 
•Initial gas phase velocities beyond this distance were 
ramped to the bulk velocity of 1.0 m/s. Spray droplet data 
were specified at radial distances out to 13.0 mm in 1 mm 
increments. The total number of droplet groups measured at 
this axial location was 102 as the large droplets were not 

found near the center line. 

According to McDonell, the error in vapor concentra- 
tion measurement was about 10 %. The error in droplet 
measurements is higher, especially in the initially dense 
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spray. The sum of the measured methanol vapor and droplet 
fluxes only accounted for 23 % of the known methanol flow at 
the first axial measurement location. The fraction of the 
measured methanol flux goes up to 85 % at the last measure- 
ment location. The PDI measures the number of fringes 
crossed and does not accurately measure flows in very dense 
sprays. Droplets are going through the measurement zone at 
various angles, velocities, shapes, and sizes. The measure- 
ments are invalidated if there are multiple droplets detect- 
ed in the measurement zone. As a rough approximation, the 
corrected droplet count rate is usually taken as the sam- 
pling rate. The error in measuring small droplets in dense 
sprays is larger as there exists a higher probability that 
multiple droplets are in the measurement zone. The ratio of 
the sample rate to the validated rate is around 5 for the 
smallest droplets at the first axial measurement location. 
This ratio drops to 1.1 for the largest droplets in the 
outer less dense region of the spray. 

Four sets of spray droplet input data were used in the 
predictions. For the first calculation the droplet data 
rate was incremented at each radial location so that the Mc- 
Donell's flux rate at that radial location was achieved. 
Then, an overall final correction was done for the flux 
rates so that the total methanol flow rate was correct. 

This droplet rate correction seemed to bias the numbers 
towards the larger droplets. Larger droplets evaporate more 
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slowly than smaller droplets as the surface to volume ratio 
is smaller for the large droplets. The calculated methanol 
vapor flux was lower than that measured by McDonell for the 
last two stations, as shown in figure 5.11. In the second 
prediction, the sampling rate was simply multiplied by an 
overall value to give the total methanol flux. This slight- 
ly increased the vaporization rate. For the third predic- 
tion, the actual probe area values were obtained from Mc- 
Donell. These probe volumes changed from location to loca- 
tion. This information allowed for a more accurate measure- 
ment of droplet flux. The calculated vaporization was simi- 
lar to that obtained in previous predictions. In the fourth 
prediction, the overall mass correction was weighed for each 
droplet size according to the ratio of the sampling rate to 
the validated measurement rate. This calculation had the 
largest number of small droplets and predicted the highest 
vaporization rate. It was hoped in using different initial 
spray distributions that the velocity field predictions 
would improve. As shown by figures 5.12-15 this did not 
happen. The use of different droplet data did not drasti- 
cally change the axial velocity profiles. Thus, the error 
in predictions does not appear to be caused by the initial 
spray distribution, or the spray model. Some of the error 
may be due to the initial radial gas phase velocities. 
McDonell 's data for these velocities is sparse. Near the 
atomizer, measurements were taken only very 4 mm 
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Figure 5.12 Axial Velocity at 15 mm 
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distance of 12 mm, the reported radial velocity increases to 
over 3.2 m/s and then decreases to -1.547 m/ s . The turbu- 
lence is extremely high. Axial velocities vary as much as 
5% between the axial-radial and the axial -azimuthal velocity 
sweep measurements. In the predictions, the initial high 
radial velocity components caused a rapid decay in radial 
and axial velocity components. To check the radial velocity 
components, another calculation was done with much smaller 
radial velocity component near the atomizer. This produced 
axial velocity profiles with higher peaks. Also, the peaks 
occurred along the centerline, as is shown by most of the 
measured velocity profiles. With the reduction in radial 
velocity, the predicted velocity profile almost perfectly 
matches the measured velocity profile at 50 mm, but, at 100 
mm the predicted profile is flatter. The predictions appear 
to be suffering from too much diffusion. Various curvature 
effects can be applied to coefficients in the turbulence 
model, but, with the rather small amount of curvature in 
this flow, these corrections would do relatively little to 
improve predictions. 

Predicted methanol vapor concentration is shown in 
figure 5.16 for an axial distance of 15 mm. The contours 
are remarkably similar for all of the predictions at this 
axial location. Experimental measurements were not given at 
this location. A comparison of experimental and predicted 
contours at 25, 50 and 100 mm is given in figures 5.17-19. 
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Figure 5.18 Methanol concentration at 50 mm 
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All of the predicted profiles are flatter and wider than the 
measured profiles, with the highest peak shown by droplet 
distribution number 4 which has the largest number of small 
droplets. The predicted methanol contours are diffusive, 
suggesting a deficiency in the numerical differencing model. 
Correcting this deficiency is beyond the scope of this work. 
The outer portions of all of the predicted contours are 
similar. McDonell states that some of the methanol concen- 
tration measurements are at a saturation state along the 
center-line. Predicted methanol concentration initially 
drops along the centerline. The maximum calculated concen- 
tration is never as high as the peak experimental measure- 
ments. This suggests that the spray model may be somewhat 
deficient in regards to the saturation concentration. An 
off-center peak is shown by all of the predicted contours at 
100 mm while the measurements show off-center peaks at 50 
and 100 mm. In the predictions, an axial location is even- 
tually predicted in the flowfield where all droplet groups 
are continually spreading outward. Even if a droplet group 
is started with a negative radial velocity, the calculated 
droplet trajectory will eventually cross the centerline, be 
reflected in the calculation and become positive. Thus, the 
calculation should predict an off-center methanol concentra- 
tion peak as the evaporating methanol droplets move further 
from the center-line. Some predicted and measured drop- 

let axial velocities at axial distances of 15, 25, 50, and 
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75 nun are shown in figures 5.20-23, respectively, for the 
third prediction (probe area taken into account, no individ- 
ual weighing of error measurement) . The data for the other 
predictions should be almost identical, as the largest 
difference should be the gas-phase velocities which was 
shown to be almost identical for all the changes in spray 
input. The predicted velocities are shown as curves while 
the measurements are shown as symbols. The predictions were 
done with initial droplet sizes of 7.4, 15.5, 25.5, 35.5, 
45.5, 55.5, 68.0, 83.0, 98.0, and 120.0 \m in diameter to 
correspond to the ten size ranges used in the measurements. 
The profiles at 15 mm are very good and remarkably smooth. 
This tends to verify the initial conditions used for the 
calculation. The predictions for the large droplets at the 
next axial measurement location are also very good. The 
prediction for the 31 to 40 jun droplet groups show some 
oscillation. Some of the droplet groups are crossing each 
other. One group with lower axial velocity is crossing the 
path of a droplet group with higher axial velocity, result- 
ing in a discontinuous profile when plotting the velocities 
of the droplet groups. The further downstream, the more 
opportunity for the predicted droplet groups to cross paths. 
However, all of the droplet groups do approach the gas phase 
velocity as a function of their size and speed. This tends 
to smooth the velocity profiles, canceling the first effect. 
Also, this effect should not have much of an effect on the 
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Figure 5.23 Droplet and gas-ph.ase axial velocity at 75 mm 
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concentration prediction, which is more important than the 
prediction of each droplet group's velocity. At an axial 
distance of 50 mm the predicted velocities for the two 
smallest groups shown are nearly the same as the gas-phase 
velocities while the measured profile shows the same behav- 
ior at less than 28 mm radius but shows a second outer peak 
with nearly the same velocities as measured for the 31-40 
and 51-60 SMD droplet group measurements. In checking the 
predictions more closely, it was found that no droplets in 
this size range were found beyond a radius of 18 . 5 mm at 
this axial location. The initial measurements at 7 . 5 mm do 
not show a secondary peak. There is' no gas-phase fluid 
dynamic structure or effect that would cause this. The 
only possibility is droplets from another droplet group. In 
the predictions, only the smallest droplets significantly 
evaporate, while the larger droplets only lose a small 
fraction of their size. The smallest droplets eventually are 
treated as totally evaporated when they reach a certain 
size. This is done because the fixed droplet time step 
could produce non-physical results. The use of smaller 
droplet time-steps would allow tracking the droplets longer 
but this would increase the number of computations and 
memory needed dramatically. The predictions do not show any 
droplets losing enough mass to be included in a lower size 
droplet group within the first 100 mm. The predictions show 
that some of the original 15.5 |im droplets are reduced to 
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10.5 |im in size when they pass the last measurement location 
of 100 mm. Larger diameter droplet groups are reduced in 
diameter even less. The largest droplets only lose about 2 
(im in diameter in 100 mm. The predictions are started using 
10 different droplet sizes. Each inlet droplet represents 
an average of a very large number of droplets spanning at 
least 10 |im in diameter. Conceptually, a group of droplets 
near the smaller end of a droplet bin size range would drop 
into the next size measurement bin at a downstream location. 
The secondary peaks measured later in the flowfield are most 
likely part of a group of droplets having larger diameter 
and velocity near the injector. This effect could be dupli- 
cated in a calculation by increasing the number of initial 
sizes. The smaller size would evaporate enough to drop into 
the next size range during the simulation, providing data 
for the smaller-diameter group at that location. However, 
the prediction of these velocities is not a requirement for 
accurate calculation of methanol concentration. The pre- 
dicted droplet is still there, producing vapor, and other 
effects. It's just that the spray model works with a dis- 
crete size and the measurements deal with a range of sizes. 
The production of methanol vapor is largely unaffected if 
additional droplet groups are added with the appropriate 
averaging for the diameters. More droplet groups would help 
along the calculation axis of symmetry where calculated 
droplet trajectories eventually all diverge. The spray 
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model interpolates between droplet sources for individual 
grid cell source terms. With no or few source terms near 
the axis of symmetry, the predicted methanol concentration 
can be depressed. 

There is another difference between predictions and 
measurements. The predictions show large droplets up to and 
beyond a radial distance of 60 mm at an axial distance of 50 
mm. The measurements were stopped at an radial distance of 
48 mm. The measurements at this location show a fairly high 
droplet rate. From these numbers, it may be deduced that 
there are many droplets beyond 48 mm. This same defect is 
also seen at an axial distance of 75 mm, where the calcula- 
tion shows droplets up to 80 mm radius and the measurements 
are cut off at 60 mm. This defect shows up also in bar 
graphs of the droplet flux. 

In the measurements, the number of droplets crossing 
the probe area are counted. A radial profile of these 
numbers is somewhat Gaussian. In the predictions, the 
droplet flux must be converted into a number of droplets 
entering the calculation domain in a certain time-step. The 
number of droplets crossing an area must be converted into a 
set of droplet groups with specified numbers, sizes and 
velocities. The Lagrangian predicted radial location of 
droplet groups changes as the groups move downstream. The 
predicted radial location does not vary uniformly as does 
the experimental measurements which are taken at specific 
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locations. Also, in the predictions, the number of droplets 
in a group does not change, unless that group evaporates (is 
taken out of the calculation). Thus, the measurements of 
density and the predictions of droplet location cannot be 
directly compared. In this case, the predicted data will be 
processed to resemble a normalized droplet density to com- 
pare to the measured droplet fluxes. This was done by 
dividing by a pseudo-area for the droplet group and then 
non-dimensionalizing this by the total number of droplets in 
that size range. A minimum area was specified so that a 
unrealistically large number would be calculated at small 
radii. Also, if more than one group was calculated at 
nearly the same radius, the groups were sometimes combined 
in the plots. The measured data are also normalized. That 
is, the sum of all normalized droplet densities adds up to 
1.0. Accordingly, the data obtained are a somewhat graphi- 
cal representation of the calculated position of the droplet 
groups and the fraction in each group. The droplet distri- 
butions at 15, 25, and 75 mm are shown in figures 5.24-31. 
Both the measured and predicted droplet densities are some- 
what Gaussian and show spreading. The numerical predictions 
generally keep the same shape through the calculation until 
droplet trajectories cross each other. At this point, the 
predicted bar graphs become somewhat ambiguous. Droplet 
group fractions should probably be combined in some of 
plots. In the predictions, the effect of an increase in 
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Figure 5.24 Measured droplet density data 
for 11 - 20 micro-meter droplets 



Figure 5.25 Calculated Droplet density for 
10.5 - 20.5 micro-meter droplets 
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droplet density is seen as a decrease in distance between 
droplet groups, not as an increase in the number of drop- 
lets. An increase in droplet density is intuitively contra- 
dictory for this flowfield, and numerical predictions of 
such effects most likely represent slightly incorrect ini- 
tial velocities. However, the somewhat inaccurate calcu- 
lation of numerous droplet densities does not have a large 
sf fsct on the prediction of methanol vapor concentration. 

It would be computationally impossible to model each drop- 
let. The numerical predictions do not show as much spread- 
ing as do the measured droplet densities, except for the 
largest diameter droplets. The spray prediction model does 
not take into account interaction between turbulent eddies 
and droplet groups. Droplet turbulent dispersion models 
take this effect into account. In this model, droplet 
interactions with the turbulent eddies are modeled by vary- 
ing the instantaneous gas phase velocity surrounding the 
droplet by using statistical fractions of the turbulent 
fluctuation velocity calculated from local turbulence quan- 
tities. This would tend to increase the droplet spreading. 
However, to do this the number of droplet groups must be in- 
creased by at an order of magnitude in order to maintain 
meaningful results. This would increase calculation time 
also by an order of magnitude. This is impractical for this 
work as testing turbulent combustion models will also sub- 
stantially increase calculation time. 
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The largest diameter droplet density graphs show simi- 
lar spreading for both experimental and numerical predic- 
tions. The larger droplets are affected less by gas-phase 
drag, having a lower surface to volume ratio. Thus, the 
larger droplets tend to spread further. The full extent of 
large droplet spreading was not found since the measurements 
were terminated too early. 

The current droplet spray model produces good results 
for methanol concentration. Droplet velocities can be sur- 
prising well predicted. There is an error in the amount of 
droplet spreading. However, there isn't a large enough 
variance from experimental measurements to justify a turbu- 
lent dispersion spray model, with its substantially increa- 
sed computational costs. 
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5.2.3 - Spray Combustion Calculations 

Both spray modeling and spray measurement capabilities 
have markedly improved in the last 15 years. In a 1980 
paper, El Banhawy and Whitelaw 113 assumed two spray distri- 
butions. Also, the droplet vaporization model used in the 
calculations did not take into account the effect of combus- 
tion. The capability now exists to measure droplet spray 
distributions for size, number flux, and velocity. Spray 
combustion data is being assembled from simplified research- 
type combustors. Gaseous data is also needed to validate or 
benchmark spray combustion models, as the gas phase must be 
accurately predicted in conjunction with the liquid phase. 
Thus, measurements of gaseous velocity, species concentra- 
tion and temperature must also be taken. There is some 
progress at measuring species and temperature using non- 
intrusive optical methods, but, this is a developing and 
specialized discipline. Measurements of temperature and 
•species can have substantial error. Spray measurements at 
the inlet of the combustor have the greatest impact on 
numerical predictions. These measurements are also the most 
difficult . The flow is optically dense close to the spray 
nozzle. Thus, the more accurate spray measurements are at 
substantial distances from the spray nozzle. The accuracy 
of inlet conditions for spray combustion leaves room for im- 
provement. The gaseous species concentrations may have to 
be inferred from the difference of liquid fuel delivered to 
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the combustor and measurements of fuel spray flux. 

A partial database for spray model validation is given 
in Edwards and Rudoff 114 . An optically dense zone of rough- 
ly 25 mm diameter and 25 mm long is shown in one of the 
figures. Spray data are given at an axial station of 50 mm. 
Velocity vector diagrams for the gas phase and spray show 
that the reverse flow along the combustor centerline reaches 
no closer than 25 mm from the spray nozzle. Droplet data at 
the first two measurement locations are marked as biased. A 
central recirculation zone is caused by swirling inlet flow. 
The flow is not contained. Maximum velocities within the 
burning flow are of the order of 15 m/s. The rest of the 
flow appears quiescent. This type of open, recirculating 
flow is particularly difficult to accurately predict. Phase 
Doppler Particle Analyzer or PDPA is used in the measure- 
ments. In a 1988 paper 115 , the PDPA measurement method was 
compared to another method. It was found to agree within 
15%. In a 1993 paper 116 , a correction scheme is demon- 
strated. It was concluded that this correction scheme was 
difficult to employ in swirling 3-D flows. A new correction 
scheme is being worked on. 

In another 1993 paper, Presser et. al. 117 present some 
non-burning and burning data for swirl numbers of 0.0, and 
0.53. Laser sheet pictures are shown for swirl numbers of 
0.0, 0.39, 0.53, and 0.76. The gas phase flow velocities of 
the open combustor are low, giving a moderate sized central 
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recirculation zone. The maximum axial velocity is about 3 
m/s. Maximum axial spray velocity is over 15 m/s. Droplet 
axial velocities seem to be all positive for a swirl number 
of 0.53. Increasing swirl increases the radial displacement 
of the spray. Temperature and species concentration mea- 
surements are not given. 

In another 1993 paper, Ghaffarpour and Chehroudi 118 
present temperature and droplet data for a hollow cone 
kerosene spray. Droplet characteristics were reported for 
mean velocity and mean droplet diameter. Negative mean 
axial droplet velocities are shown along the combustor 
centerline at distances of 40, 55, arid 70 mm for the non- 
burning case. The first measurement location was 15 mm from 
the injector. Mean axial droplet velocities at this station 
were all positive, although the velocities were lower for 
the burning case. The burning case shows considerably more 
radial spreading of the fuel spray. The only gas phase 
velocities reported were from the non-burning case, which 
indicated a central recirculation region. Measurements of 
burning and non-burning droplet flux are nearly identical, 
indicating negligible droplet vaporization and/or burning. 
Conversely, temperatures within the hollow cone area were 
fairly high, about 500 degrees C at 15 mm and 600 degrees at 
25 mm. The authors infer that this is due to a downward 
movement of the central recirculation zone for the burning 
case, but this flow needs to be properly seeded to measure 
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velocities in this region. Before this case can be used for 
spray model development, the gas phase velocities must be 
accurately measured and droplet statistics should be speci 
fied for various droplet size groups, not just for mean 
diameters. Using mean quantities limits the number of drop- 
let groups at the 15 mm location to about sixteen. 

The combustion case of Bulzan [108] was used to bench- 
mark the combustion spray model. The flow configuration is 
identical to the isothermal case studied earlier, except 
that the injector was supplied with liquid heptane fuel. 
Temperature, gaseous and liquid velocities were measured. 
Data were meticulously taken and repeated to minimize exper- 
imental error. Great pain was taken to eliminate 3-Dimen- 
sional effects, such as locating the swirler for the co- 
flowing air-stream 140 mm upstream of the injector face. 

The combustion experiment does not have combustor walls. 

This makes it simpler to get optical measurements, but makes 
it difficult for the combustion modeler, as the flow en- 
trains ambient air. Runs were made with entraining boundary 
conditions, but this modeler found it best to specify bound- 
ary conditions at measured radial locations. The average 
radial velocities are rather low but show large variation. 
This modeler was not able to duplicate such high turbulence 
with low velocity. Initial calculations were performed 
using 70 axial gridpoints and 55 gridpoints in the radial 
direction. Measurements taken 2 . 5 mm from the face of the 


spray injector were used as initial conditions. These data 
are much closer to the combustor spray nozzle than reported 
by others. The droplet groups are differentiated by droplet 
size ranges and location. Close to the injector, droplet 
measurements were taken at 1 mm intervals out to a radius of 
about 13 mm, where the data rate became negligible. The 
droplet data were separated into 10 size ranges or bins. 
These bins are then characterized by average diameters. The 
average group diameter range is from 6.893 to 122.65 micro- 
meters. These data were reduced to 79 droplet groups as 
some velocity data were unreliable due to the scarcity of 
data, or the amount of fuel spray in the group was judged as 
insignificant. Obviously, the computer time needed to 
calculate the spray and its influence on the gas phase will 
roughly scale with the number of groups. Data were taken as 
a number flux at specified radius. For calculation purpos- 
es, the flux data were multiplied by an appropriate area to 
get the total number of drops in a group. The total spray 
mass was calculated and this was significantly different 
than the measured fuel flow rate to the injector, which was 
originally taken as 0.38 g/sec. The difference was input as 
a combination of fuel and combustion products. Originally, 
it was planned to take probe samples for discerning gaseous 
species concentrations. Bulzan now feels that sample probes 
would distort the flowfield producing inaccurate results. 
Measured values of temperature are low in the dense spray 
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due to droplet impingement on the probe. Gas sample probes 
would also be greatly affected by the spray. The only data 
with which to infer gaseous mass fractions are the total 
rate of fuel flow and temperature measurements. The temper- 
atures in this area were inferred from neighboring tempera- 
tures not exhibiting droplet cooling. To come up with 
appropriate inlet mass fractions, a sub-program was used to 
calculate temperature as a function of increasing fuel mass 
fraction. A fixed final ratio of unburned to burned fuel 
was used while keeping the total enthalpy constant. The 
numbers were iterated/modified upon to get the total overall 
fuel rate. As initial rough grid calculations produced 
higher than measured temperatures, the heating value of the 
fuel was dropped by 10 % for better agreement. Initial 
calculations were done assuming no correction to droplet 
number flux. In McDonnel's thesis, it was concluded that 
most of the droplets were not measured in the area of dense 
spray close to the injector. If multiple drops are detected 
within the measurement zone, the data are discarded. The 
measurement of mass fractions were determined to be more 
accurate than that of droplet flux, so the droplet flux was 
corrected. The droplet flux closest to the injector was off 
by a factor of almost 4. Bulzan didn't report species 
concentrations which would have enabled correcting inlet 
droplet number flux. Also, the data are further complicated 
by combustion, widely varying temperature, and very rapid 
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evaporation. Several 70 by 55 grid calculations were made 
with different spray flux corrections. It was found that 
the best agreement was found by assuming that about 35 % of 
the fuel and fuel in combustion products were gaseous. 
Temperature was significantly overpredicted. Since the 
rougher grid calculations were made, it was found that the 
total fuel flow rate was actually 20 % lower. Succeeding 
calculations used a finer grid and a corrected overall fuel 
flow rate. The finer grid essentially doubled the number of 
radial gridpoints in the fuel injector zone. A velocity 
vector diagram of a portion of the flowfield is shown in 
figure 5 . 32 . The flowfield structure is simpler than for 
the isothermal case. In the combusting case, the recircula- 
tion bubble is closed within the confines of experimental 
measurement limits. Axial velocities are all positive 200 
mm from the combustor injector. In the isothermal case, 
reverse flow existed out beyond the last measurement sta- 
tion, and modeling proved to be rather difficult. In addi- 
tion, the small toroidal recirculation bubble seen near the 
fuel injector and co- flowing air stream in the isothermal 
case isn't measured or predicted for the combusting flow 
case. Calculation of entrainment is still a problem for the 
combusting case. 

Measured centerline axial velocity next to the fuel 
injector is negative. This reverse flow is part of a large 
central recirculation bubble caused by high inlet swirl. 
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Theoretically, the reverse flow should end very close to the 
fuel injector. The closest Dr. Bulzan was able to perforin 
measurements was 1.0 mm from the injector. At this loca- 
tion, the axial velocity was non-negative, but the spray 
measurement data were in question. It may be possible to 
perform calculations starting from the injector face, but 
the fuel stream breakup would then have to be modeled, which 
cannot be done at the present time. In this work it was 
desired to see how good the spray model was, without cou- 
pling it to additional modeling to develop or back calculate 
suitable droplet distributions. It was decided to start 
with the best possible measured droplet distribution and see 
how well the combustion spray model predicted successive 
axial spray distributions. 

The numerical simulations started 2.5 mm from the 
injector face. Plots of the experimental and numerically 
predicted axial velocities are shown in figure 5.33-38 
corresponding to axial distances of 5.0, 10.0, 20.0, 50.0, 
100.0, and 200 mm from the fuel injector. In taking the 
experimental measurements a complete traverse was taken 
across the combustor. The measurements on one side of the 
combustor are plotted with negative radii. Numerical pre- 
dictions were 2-D axisymmetric , thus, the numerical results 
are simply reflected for comparison purposes. The axial 
velocity predictions are quite good beyond the region of 
fuel injector flow. Surprisingly, there is little 
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difference between numerical predictions despite the model- 
ing changes . 

The inlet reverse flow velocity specified next to the 
centerline was about -23.0 m/s. The predicted centerline 
velocity at an axial distance of 5 mm is about 10 m/s, while 
the measured centerline velocity is about -30.0 m/ s . Exper- 
imental centerline axial velocity increases while the simu- 
lations predict a decrease in centerline velocity. This 
poor prediction pattern is also seen for peak axial velocit- 
y. The numerical simulations predict peak axial velocity of 
approximately 10 m/s at an axial distance of 5 mm, while 
measured peak axial velocity is almost identical to the 
inlet peak velocity of 35 m/s. The numerical simulations 
are poorly predicting a critical zone of the combustor. 

This area of the flow contains the highest concentration of 
fuel spray. The gaseous velocities should be correctly 
predicted in order that spray predictions are also correctly 
-predicted. Many additional calculations were done modifying 
turbulent inlet conditions and modifying the turbulence 
model in an attempt to better predict reverse flow along the 
combustor centerline and peak axial velocity. These simula- 
tions did not appreciably improve the prediction of axial 
velocity close to the injector, but mostly shifted the axial 
velocity peak in the radial direction. 

Experimental and predicted peak axial flow velocities 
at an axial distance of 10 mm are smaller than at the 5.0 mm 
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station. Predicted peak axial velocities are roughly one- 
third of the measured velocities. The simulations predict 
almost complete merging of the air-assist fuel injector flow 
with the surrounding co- flowing air-stream. The simulations 
predict a larger radial displacement of the fuel-air injec- 
tor flow than measured. Predicted reverse flow centerline 
velocities are less than half of measured velocities. 

The experimental data taken at an axial distance of 20 
mm shows distinct peaks due to the fuel-air injector flow. 
The simulations are nearly identical except in the central 
recirculation zone. The simulations predict a single veloc- 
ity peak. Qualitatively, predicted axial velocities are 
fairly good, except in the central reverse flow region. 
Predicted reverse flow velocities are about half of the 
measured velocities. 

At 50 mm, the axial velocity predictions seem to im- 
prove. Single velocity peaks are shown for predictions and 
measurements. The simulations slightly underpredict peak 
velocity. The quality of the axial velocity predictions in 
the central recirculation zone are much improved from previ- 
ous axial stations. The finest grid shows the lowest cen- 
terline velocity. 

Peak velocity is well predicted at an axial distance of 
100 mm. There is little difference between the two numeri- 
cal simulations. The width of the predicted axial velocity 
profiles are noticeably thinner than the measured profiles. 
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This implies the simulations may be incorrectly predicting 
the total mass flow in the axial direction. This is proba- 
bly due to the complexity of calculating an open recirculat- 
ing flow which is entraining ambient fluid. 

At 200 mm, the simulations predict nearly identical 
profiles. Centerline axial velocities are all positive. 

The recirculation zone closed somewhere between 100 and 200 
mm downstream of the fuel spray nozzle. Axial velocity is 
well predicted from a radius of 5 cm to 15 cm. The simula- 
tions predict parabolic -type axial velocity profiles that 
are symmetric about the peaks. The experimental profile is 
noticeably non- symmetric around the peak. Large axial 
velocity gradients are measured about the combustor center- 
line. The simulations predict a much lower gradient. 
Qualitatively, the large velocity deficit at the combustor 
centerline is well predicted. The velocity deficit along 
the centerline will disappear much more quickly in the 
experimental profile than in the predicted profiles. The 
local minimum experimental axial velocity isn't located 
exactly at the combustor centerline. The experimental axial 
velocity profile is slightly shifted to the right. This is 
probably due to a slight asymmetry in the flow. 

The radial velocities are shown in figures 5.39-44. 

The peak radial velocity is underpredicted by both simula- 
tions. However, the radial velocity peak is better predict- 
ed than the axial velocity peak. The fuel-air flow 
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Figure 5.39 Radial Velocity at 5.0 mm. 



Figure 5.40 Radial Velocity at 10.0 mm. 
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expansion due to combustion seems to occur more in the 
radial direction. The simulations predict a larger radial 
expansion of the fuel-air injector flow than measured. The 
predicted width of the fuel-air injector flow is nearly 
twice as wide as measured. The predicted radial over-expan- 
sion destroys the large radial velocity peak at successive 
axial locations. At 10.0 mm, the predicted radial velocity 
peaks are about twice as far from the centerline as the 
experimental peaks. The predicted peak radial velocity is 
less than one-fifth of the experimental peaks. Radial 
velocities beyond a radius of 4 cm are well predicted. 

At 20.0 mm, both the numerical simulations predict a 
single radial velocity peak. The flow from the fuel-air 
injector has completely merged with the co- flowing stream in 
the simulations at this axial location. The experimental 
data shows separate velocity peaks. The experimental veloc- 
ity peaks are about four times higher than predicted peaks. 
Radial velocity in the co- flowing stream is well predicted. 

The experimental data shows a single radial velocity 
peak at the 50.0 mm station. Flow originating from the 
fuel-air injector and the co- flowing air-stream have com- 
bined. Qualitatively the simulations correctly predict the 
radial velocity at this station. The simulations under- 
predict the velocity peak by about 50 percent. The radial 
velocity in the central recirculation zone is off by at 
least a factor of four. The negative velocities near the 
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combustor centerline indicate the recirculation zone is 
starting to close. 

The radial velocity profiles for the entraining flow 
and outer co-flow are well predicted for the 5.0, 10.0, and 
20.0 mm stations. However, the portion of the outer profile 
that is in agreement, grows progressively smaller with axial 
distance. At 50.0 mm, the agreement in the outer flow is 
the smallest. The width of the predicted profile is under- 
predicted. There is the thought that the initial radial 
velocities may be over-specified, causing an over-rapid 
broadening of the injector flow. This researcher did a 
calculation specifying reduced radial velocities for the 
injector flow. This modification remarkably improved the 
calculation of the peak velocity flow from the injector. 
However, there is little evidence to support such a change 
in radial velocities, other than it works. Another possi- 
bility of error is the flow direction to the computational 
grid. The injector flow is roughly 45 degrees to the grid 
lines, which is known to cause accuracy problems. It would 
be interesting to implement an improved or higher order 
numerical differencing scheme in the model. 

At 100 mm, the numerical radial velocity predictions 
are all negative. That is, there is inflow across the whole 
combustor, no expansion. The experimental data shows some 
expansion at larger radii. At the 200 location, the pre- 
dicted radial velocities are very close to zero. 
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sured velocities at this station are generally larger. 

There is significant non-symmetry for the 100.0 and 200.0 
experimental data near the centerline. The velocity should 
go to zero near the centerline but doesn't. Part of the 
problem is the sign of the radial velocity. It is believed 
that the flow is slightly off to one side, which causes 
problems in interpreting the sign of the radial velocity 
component. The experimental axial velocity profile exhib- 
ited sharp gradients near the centerline for the 100.0 and 

200.0 stations, thus, there should be significant radial 
velocity in this area. 

The plots of tangential velocity are shown in figure 
5.45-50. As in previous figures, the part of the flow 
agreeing least with the experimental data is the injector 
flow. The numerical predictions at 5.0, 10.0, and 20 mm 
show a large drop in tangential velocity for the injector 
flow. The inner profile is worst predicted at the 10.0 mm 
station, where it is a quarter of what it should be. The 
predicted inner profile actually gains some strength at the 

20.0 mm station. The experimental data for the injector flow 
exhibits quite a bit of variability and asymmetry. Between 
the 5.0 and 10.0 station, the experimental velocity peak 
changes sides! The difference in injector flow velocities 
from one side to the other is of the order of 25%. 

The outer velocity profiles are well predicted at the 
5.0, 10.0, and 20.0 locations. Then the width of the 
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profile is off, but predictions significantly improve! At 
50.0, 100.0, and 200.0 mm the prediction of tangential 
velocity is good. The velocity peaks are only slightly 
under-predicted. The width of the profile are underpredict- 
ed, caused by a sharper drop in tangential velocity than is 
measured. The inner profile is underpredicted at the 100.0 
and 200.0 locations. This was also seen for the axial 
velocity profile. The experimental data indicates stronger 
radial inflow near the centerline. Viscous effects are 
acting on the on the strong velocity gradient near the 
combustor centerline. The measured turbulence is quite high 
at this location. At 350.0 mm, the experimental axial 
velocity profile only shows a slight deficit along the 
centerline. 

Plots of temperature profiles are shown in figure 5.51- 
56. The agreement on the left side of figure 5.51 is good. 
The left side agreement is noticeably better for the 5.0 
•through 20.0 mm stations. It appears that initial combus- 
tion is better on the left side for these stations. This is 
the side that the flow seems to shift or expand away from. 
The predicted profiles have rather sharp ridges at the outer 
edges of the high temperature region that are probably non- 
physical . The depressed temperatures due to fuel droplet 
impingement on the thermocouple probe are evident . The 
model predicts lower temperature in the fuel spray zone. 

Both experimental and predicted temperature profiles show 
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very large gradients at the edges where combustion is not 
occurring. The finer-grid calculation predicts higher 
temperatures than the coarse grid. At 50.0 mm, the experi- 
mental temperature profile is fairly flat across the combus- 
tion zone. The coarser grid exhibits two small humps, while 
the fine grid has very prominent peaks. The peaks are much 
smaller at the 100 mm station. The width of the temperature 
profile is slightly underpredicted at the 100.0 mm station. 
The gradient at the edge of the high temperature profile is 
overpredicted. This is probably due to the underprediction 
in radial velocity. At 200 mm, the experimental and pre- 
dicted profiles show a rather bell shaped curve. The exper- 
imental profile shows a little more spreading, again proba- 
bly due to the radial velocity. 

Fuel Droplet Comparisons 

The droplet axial velocity comparisons are shown in 
figures 5.57-60. Droplets were not measured past the 50 mm 
station. The experimental data are shown for negative radii 
in the figure, and the fine grid numerical predictions are 
shown for positive radii. The solid line is the gas phase 
velocity shown for comparison purposes, as the droplets will 
eventually relax to the gas phase velocity due to drag. The 
experimental gas phase velocity plots were not smoothed. 

The different symbols denote the different size droplet 
groups. Only 6 of the 10 droplet sizes are plotted. At 5 . 0 
mm, the numerical calculation of spray droplet axial 
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velocity is slightly underpredicted. This is caused by the 
under prediction of the gas phase axial velocity near the 
vicinity of the fuel injector. The larger difference be- 
tween the predicted gas phase and droplet velocities results 
in higher drag forces, and a quicker reduction in droplet 
velocities. The smaller droplets generally track the gas 
phase velocity better. The shape of the droplet velocity 
profile is good, except, the droplet groups show greater 
radial spreading than experimentally measured. The experi- 
mental data at 5.0 mm indicates many particles near the 
combustor centerline. Many of these particles show negative 
axial velocities! The predictions do not show the same 
patterns. There are few predicted droplet groups near the 
centerline, and only one with negative axial velocity. 
Gaseous fuel is provided by evaporating droplets. The 
combustion of this fuel results in expansion and possible 
increase in reverse flow velocity. The experimental data 
shows an increase in negative axial velocity at this loca- 
tion, while the predictions show a decrease. The spray 
model is not picking up on possible phenomena to explain 
this effect. In setting up the droplet group input data, 
droplet groups with negative axial velocities were ignored. 
It would be useless to input this data. The spray flowing 
towards the injector significantly contributes to gaseous 
fuel mass fractions, but without doing some species measure 


merits, it is unknown. 


At 10.0 mm, The droplet axial velocities have de- 
creased further. The predicted gas phase velocity does not 
even show a peak corresponding to the air-assist stream from 
the nozzle. The difference in predicted velocity for the 
two phases is much higher than the experimentally measured 
velocity difference. This will lead to further loss in 
droplet axial velocity. The predicted droplet velocities 
exhibit a much sharper peak than measured and show greater 
radial spreading. The experimental droplet data shows 
grouping of the data according to droplet size. Droplet 
groups having the same average diameter seem to plot along 
smooth curves. This isn't as obvious in the numerical 
predictions, especially, later in the flowfield, due to 
scarcity of particles. 

The number of predicted droplet groups decreases as 
the droplets are burned. There are only 5 droplet groups 
plotted at the 20.0 location, and two at the 50.0 mm sta- 
tion. This is far too few show droplet size grouping. It 
is also questionable if this small number could cause any 
kind of local gas phase velocity peak. The experimental 
data shows good correlation between the injector gas veloci- 
ty peak and the placement of the droplets. The simulations 
don't adequately predict the air-assist nozzle flow. To 
help validate the spray model, it would be better to have 
significantly more droplet groups. Unfortunately, the 
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number of droplet groups was rather limited to begin with! 
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The current spray model significantly taxed computer time 
resources. The needed computer time would scale directly 
with the number of droplet groups. The experimentalist can 
widen the measuring probe volume to pick up more droplets 
where the spray is less dense. There isn't a numerical 
modeling counterpart to this. If one mentally averages the 
droplet axial velocities, the simulation seems about at 
third of the experimental average for both the 20.0 and 50.0 
measurement stations. Intuitively, this should have a 
severe effect on the spray vaporization. Dr. Bulzan inte- 
grated his spray flux measurements across the flowfield. 
Correcting these measurements gives 21.3 % of the metered 
flow rate at 2.5 mm, 26.8 % at 5, 50.0 at 10, 23.1 at 20, 
and 2.6% at 50.0 mm downstream of the nozzle. It is obvious 
that much of the droplet data are not picked up at the 2.5 
and 5.0 stations. The predictions presented here used 65% 
of the metered fuel flow for droplet input at 2.5 mm. The 
•simulation predicted 51.7 % of the metered flow in the form 
of spray at 5, 24.0 % at 10, 6.3 % at 20, 0.13 % at 50.0 mm. 
Thus, the simulation does predict over-rapid fuel spray 
vaporization. However, as velocity predictions near the 
injector zone are overly diffusive, it would be inappropri- 
ate to blame the spray model. The droplet velocity is about 
0.8 of what it should be for the 5.0 and 10.0 location. At 
the 20.0 location it is roughly a third and at the 50.0 
location it is about a quarter. Assuming the amount of 


droplet vaporization is directly proportional to time, 
applying a correction for the time by factoring in how much 

I 

the gas phase velocity differs from the data and compounding 
this correction at successive axial locations to the above 
percentages gives 57.4 % at 5, 33.3 % at 10, 15.5 at 20.0 
and 1.1 % at the 50.0 location. These amounts are still too 
low. If one looks at the 20.0 location, this suggests that 
the spray is off by a third. This would then imply that all 
the nearly all the fuel should be in liquid form at the 
calculation start. This does not agree with the data, as 
inlet temperature is very high. 

The droplet radial velocities are shown in figure 5.61- 
64. At the 5.0 mm station, there is surprising little 
velocity difference between the gas phase and spray drop- 
lets. This was also exhibited by the plots of axial veloci- 
ty. The predicted gas phase radial velocity peak is only 
slightly less than the measured one. The axial velocity 
peak was a factor of three off. However, as previously 
noted, the profile is twice as wide as it should be, which 
causes additional nozzle flow dissipation. The numerical 
simulation predicts droplets at larger radii. The cause is 
a combination of the predicted axial and radial gas phase 
velocities. The calculated radial gas phase velocity pro- 
file is much broader than the experimental one, and the pre- 
dicted axial velocities are much lower than measured. Thus, 
the predicted spray spreads out further than measured. At 
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10.0 mm the predicted gas phase radial velocity peak for the 
injector flow is much lower than the measured one. The 
experimental data at this location only shows one group with 
negative radial velocity near the centerline. The average 
diameter of this group is small, which allows it to better 
track the local gas phase velocity. As mentioned previous- 
ly, the reverse flow velocity at this point is decreasing. 

As previously discussed for droplet axial velocity, the pre- 
dicted droplet groups display more radial dispersion than 
measured. At 20.0 mm, the predicted radial velocities are 
generally much smaller than the measured. There are 5 
predicted droplet groups at this axial location. Some of the 
groups are in the high temperature flow. All of the mea- 
sured groups are within the high temperature flow. At 50.0 
mm there are only two predicted droplet groups and both of 
these are in the lower temperature flow. The measurements 
display many droplet groups at the 50.0 mm station. This 
suggests that the spray model is overpredicting evaporation. 
The experimental data at this location shows much higher 
radial velocities, both in gaseous and liquid phases. The 
spray is going to provide for gas phase source terms further 
downstream as the fuel evaporates and burns. The small 
number of predicted droplet groups means that most of the 
fuel is already gaseous and thus, already burning. From 
this information, it may be inferred that predictions should 
not noticeably improve with downstream distance. This was 
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seen in previous figures of velocity end temperature . 

Tangential droplet velocities are shown in fig. 5.65- 
68. These velocities are also under-predicted, at times, 
half the measured values. The calculated droplet velocities 
are nearly the same as the calculated gas phase velocity, 
while the experimental droplet velocities can be consider- 
ably lower than the measured gas velocities. The experimen- 
tal droplet data shows noticeable stratification by droplet 
size. The smallest diameter droplets are closest in tangen- 
tial velocity to the gas phase, while the largest diameter 
droplet groups differ the most. At the 20.0 mm station, the 
droplet groups plot nearly horizontally according to droplet 
size. The predicted droplet groups show slight pattern 
according to size at the 5.0 mm station. At later stations, 
the smaller number of predicted droplet groups does not 
allow a conclusion. The low number of predicted droplets 
suggests over-rapid vaporization and/or combustion. The 
numerical calculations do not adequately predict the fuel- 
air nozzle flow. Within a short distance, the injector flow 
has been seemingly dissipated. 

The experimental data shows higher peak radial veloci- 
ties than axial velocities for the 2.5 mm station to the 20 
mm station. D. Bulzan [108] concluded that the flow was 
expanding more in the radial direction due to the strong 
central reverse flow caused by the recirculation zone. The 
poor velocity predictions for the injector nozzle flow are 
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probably not the fault of the combustion model. The overly 
large radial dissipation of the injector flow must be inves- 
tigated. This researcher feels the most probable modeling 
improvement would be an improved numerical differencing 
scheme to adequately calculate the high velocity nozzle flow 
which is severely skewed to the current cartesian grid. 
Modification of the spray or combustion models will probably 
provide little improvement. Indeed, a preliminary calcula- 
tion was done using a pdf combustion model. While there are 
many differences in the modeling than used here, the gas 
phase velocities were very similar. The largest error was in 
calculating the injector flow, which is simply where the 
fuel spray is. If the fuel source terms to the gas phase 
are incorrectly calculated, a combustion model for the gas 
phase isn't the correct fix. 
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5.3 - Monte Carlo Pdf Model Predictions 

Four cases are predicted using the hybrid pdf method. 

The stochastic part of the coding solves for enthalpy and 
species while velocity and pressure are solved using the 
SIMPLE algorithm of Patankar. The hybrid pdf predictions are 
compared with predictions made by a more traditional combus- 
tion model which is a combination global arrhenius reaction 
rate and eddy breakup or mixing model, which was also used 
in the previous spray combustion calculations. These predic- 
tions were generally optimized to obtain good comparisons 
between measurement and predictions. 

The first case compares species concentration and tem- 
perature for an enclosed diffusion combustor. In this com- 
bustor, the fuel and air supply ducts are similar in size. 
The second case also involves a diffusion combustor, but the 
geometry is much closer to an actual gas turbine combustor. 
The fuel enters through a small conical duct which is at a 
large angle to the combustor axis. The third case involves 
premixed combustion. Premixed combustion is being studied 
as a way of lowering emissions and improving fuel efficien- 
cy. As simple combustion models may poorly predict pre- 
mixed combustion, the ability of a combustion model to accu- 
rately predict both diffusion and premixed combustion is a 
desirable feature. The fourth case involves a swirling 
hydrogen jet diffusion flame. This combustor is extremely 
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interesting as temperature pdfs were measured. 

5.3.1 Non-Swirling Diffusion Combustion 

The data of M. H. Lewis and L. D. Smoot 119 was chosen 
because of the simple geometry. Fuel and air are supplied 
by concentric tubes to a large insulated cylindrical combus- 
tion chamber. The combustion chamber is 152.4 cm long and 
20.32 cm in diameter. The fuel inlet tube has an inside 
diameter of 1.6 cm and wall thickness of 0.32 cm. The fuel 
inlet tube has a blunt end. The outside radius of the air 
annulus is 2.86 cm. This type of research combustor is 
known as a dump combustor. The large increase in cross- 
sectional area is easy to model and causes a hot recirculat- 
ing flow which provides a continuing ignition source for 
combustion. Lewis and Smoot did not measure velocities, 
although it appears that Smoot and another of his students, 
Philip Smith 120 used a CFD code to calculate flow proper- 
ties . 

The fuel is city-gas which is 88.53 % methane, 7.44 % 
ethane, 2.55 % nitrogen, 1.39 % carbon dioxide, and 0.09 % 
hydrogen. In this work, the fuel is modeled by assuming 
pure methane. The experimental inlet velocity and tempera- 
ture of the fuel was reported as 21.3 m/s and 300 degrees 
Kelvin. 20.14 m/s was used in the numerical simulation so 
that the overall fuel to air ratio was stoichiometric. The 
inlet air velocity and temperature was 34.3 m/s and 589 
degrees Kelvin. Temperature was measured at axial distances 
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of 9.5 and 39.5 cm from the combustor inlet. Species were 
measured at axial distances of 9.5, 17.5, 24.6, 32.7,47.6, 
63.2, 78.5, and 137.5 cm. The calculation modeled the 
combustor as being only 1.0 meter. long, as initial pdf 
simulations exhibited numerical difficulty due to large 
aspect ratio grid cells. The grid was stretched in the 
axial and radial directions with the finest grid near the 
inlet . 

The grid has 80 axial and 40 radial gridpoints. The 
fuel jet has 5 radial cells and the air has 11 radial cells. 
A recirculation zone is predicted. The length of the recir- 
culation zone varies between the hybrid pdf and SIMPLE 
calculations. The SIMPLE simulation predicted a reattach- 
ment length of 46.5 cm while the pdf simulation predicted 
57.8 cm. 

Smith and Smoot used an assumed shape pdf model with an 
equilibrium combustion model. This model significantly 
over-predicted carbon dioxide, especially along the center- 
line. Oxygen was under-predicted in the same area. This 
implies that burning was not occurring as rapidly as pre- 
dicted by Smith and Smoot's calculation. A finite rate 
single step combustion model is used in the simulations 
reported here. The model is: 

CH A + 2 0 2 - C0 2 + 2 H 2 0 
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W, 


w. 


w. 


Where 

r Kg-mole-i 

A = 8 . 5 x 10 13 l— “r J 

m 3 sec 

The exponents used for the fuel and oxidant concentra- 
tions are 0.2 and 1.3 respectively. These are the third set 
recommended by Westbrook and Dryer [25] for Methane. The 
first two sets recommended a negative exponent on the fuel 
which resulted in some numerical difficulty with large 
reaction rates at very low fuel concentrations* 

Comparisons of fuel, oxygen, water, and carbon dioxide 
mole fraction at an axial distance of 9.5 cm are made with 
experimental measurements in figure 5.69. There is little 
difference in the calculated species concentrations at this 
axial station. The simulations do a very good job of pre- 
dicting the fuel concentration, except for over-predicting 
the peak fuel mole fraction. Up to a radius of 2 cm, the 
oxygen concentration is well predicted by both simulations. 
Beyond a radius of 2 cm, the oxygen concentration is over- 
predicted. The pdf simulation gives slightly lower oxygen 

concentrations at large radii. 

The concentration of C0 2 was found from the concentra- 
tion of fuel, oxygen, and atom balances. Using a finite 
rate reaction model for the prediction of C0 2 resulted in 
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superfluous product concentrations at high temperatures. 

While the predicted concentration of fuel and oxygen can go 
to nearly zero, the rest of the arrhenius rate term is very 
large at high temperature. Initial computation product 
concentrations were 30 per cent above the theoretical maxi- 
mum concentration. This has serious consequences for mul- 
tiple reaction step models. 

Beyond a radius of 2 cm, product concentrations are 
underpredicted by about 50%. This is about the level oxygen 
was overpredicted. The simulations underpredict combustion 
in the outer recirculation zone. There is some disparity in 
product concentration near the combustor axis. The predic- 
tions show zero product concentration at this location. 
Measurements show significant C0 2 concentration, but near 
zero water concentration. The temperature profiles for this 
axial location are shown in figure 5.70. The temperature is 
lowest near the combustor axis. The simulations predict 
•temperatures about 200 degrees cooler along the combustor 
axis. Predicted temperatures off the combustor axis are very 
good, despite the error in predicted species concentrations. 
Temperature increases with radius up to 5 cm. Beyond a 
radius of 5 cm, temperature is fairly constant. The pdf 
simulation which shows a slight increase in temperature 
towards the combustor wall. The outer wall was modeled as 
adiabatic since the combustor wall was insulated. The 
experimental profile shows temperature slightly dropping 
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Figure 5.70 Temperature at 9.5 cm 
near the wall. 

M. Nikjooy and R. So [28] did some calculations for 
this same combustor. They used a two-step finite rate 
combustion and an equilibrium combustion model. The use of 
an equilibrium combustion model implies that infinite time 
exists for combustion or that combustion is extremely fast. 
The finite rate constant was higher than that employed in 
this work and very rapid combustion was predicted. Their 
calculations predicted peak temperatures above 2000 degrees 
close to the air and fuel inlet ducts, which is not in 
agreement with the experimental data. 

Comparisons of species at an axial distance of 17.5 cm 
for the present simulations are shown in figure 5.71. The 
SIMPLE simulation best predicts the fuel concentration. The 
pdf simulation slightly over-predicts the centerline peak 
fuel mole fraction and the fraction at a radius of 2 cm. 
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The hybrid pdf simulation better predicts oxygen mole frac- 
tion. Oxygen concentration is well predicted up to the peak 
oxygen concentration at a radius of 2 cm. Beyond the peak, 
both simulations overpredict oxygen concentration. The 
largest error is at the combustor wall. Product concen- 
trations are underpredicted by the simulations as they were 
at 9.5 cm, although the hybrid pdf simulation is improving 
compared to the SIMPLE simulation. The anomaly in product 
concentrations near the combustor axis still exists. 

The species mole fractions at 24.6 cm are shown in figure 
5.72. The peak fuel mole fraction is significantly overpre- 
dicted. At larger radii, the fuel predictions are excel- 
lent. The pdf simulation again predicts higher fuel concen- 
trations. Paradoxically, the pdf simulation also predicts 
lower oxygen concentration than the other simulation. Peak 
oxygen concentration is well predicted. Oxygen concentra- 
tion is slightly underpredicted by both simulations near the 
combustor axis. Oxygen concentration beyond a radius of 2 
cm are overpredicted, but are improving compared to previous 

comparisons. Both simulations underpredict product species 

\ 

concentration, with the exception of- water near the combus- 
tor axis . The hybrid pdf simulation better predicts oxygen 
and product concentration. Carbon dioxide is detected near 
the combustor axis, but not predicted. The H 2 0 concentra- 
tion near the combustor axis is very low for the predictions 
and measurements. Predictions and measurements of H 2 0 
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of H 2 0 concentration increase towards the combustor wall. A 
peak in C0 2 concentration is measured at a radius of 6 cm, 
while the simulations show C0 2 concentration increasing 
towards the combustor wall. 

Predictions of fuel concentration along the combustor 
axis worsen at the next axial station. At 32.7 cm, both 
simulations over-predict peak fuel concentration by thirty 
per-cent. Paradoxically, predictions for oxygen are excel- 
lent, much improved from previous axial comparisons. Under 
a radius of 4 cm, the pdf simulation only slightly under- 
predicts oxygen. Much improvement is also seen in product 
species predictions. The carbon dioxide concentration is 
slightly underpredicted. Near the combustor axis, H 2 0 
concentration is overpredicted by the hybrid pdf simulation, 
but the pdf simulation best predicts carbon dioxide. As a 
single step reaction is being employed in the combustion 
models, predicted product concentrations are related. The 
experimental data does not show this dependence between 
product concentrations near the combustor axis. Thus, evalu- 
ating only one product species at this location could lead 
to incorrect conclusions about the predictive capability of 
combustion models . 

Figure 5.74 gives the predicted and corrected tempera- 
ture profiles at 39.5 cm. The hybrid pdf simulation best 
predicts temperature near the combustor axis at this loca- 
tion. However, temperature is significantly overpredicted 
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Figure 5.74 Temperature at 39.5 cm 
beyond a radius of 2 cm. At a radius of 8 cm, the error in 
predicted temperatures is 400 degrees. Experimental temper- 
atures near the combustor wall suggest some heat transfer is 
taking place. Measurements of species concentration were 
not reported at this axial location. 

Species mole fractions at an axial distance of 47 . 6 cm 
are shown in figure 5.75. Fuel mole fraction is overpre- 
dicted by over 50%. Comparing fuel profiles for previous 
axial stations, this is the station with the largest rela 
tive error. The peak fuel concentrations have dropped by 
about a factor of two from the previous axial station. The 
fuel concentration profile is also spreading. Oxygen con- 
centration is overpredicted near the combustor axis . Pre- 
dicted oxygen concentration is lower than at the previous 
station. There is significant combustion near the combustor 
axis. Significant H 2 0 and C0 2 concentrations are measured 
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and. predicted near the combustor axis. The pdf simulation 
of product concentrations are very close to the experimental 
values. As the pdf simulation overpredicts reactant concen- 
tration near the combustor axis, this shows some error 
possibly due to the one-step reaction mechanism. By not 
allowing intermediate chemical species, product species can 
be over —predicted . The conventional combustion simulation 
underpredicts product concentration near the combustor axis 
where fuel and oxidant are overpredicted. 

Reactant concentrations at an axial distance of 63.2 
cm, shown in figure 5.76, continue to rapidly decrease. 

Both simulations under-predict oxygen near the combustor 
wall. The pdf simulation predicts higher reactant concen- 
tration and lower product concentration than the convention- 
al simulation. Predicted product concentrations are much 
more uniform than the measurements. Both simulations over- 
predict product concentrations. Product concentration pre- 
dictions are worst near the combustor axis and the combustor 
wall. Except for fuel prediction, hybrid pdf predictions 
are closer to the experimental data, especially at larger 
radii . 

At the 78.5 cm station, measurements show some oxidant 
and fuel concentration. This is despite the high tempera- 
tures that must exist. The pdf simulation overpredicts 
oxygen and fuel. The conventional simulation predicts a 
very slight amount of fuel near the combustor axis and some 
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oxygen towards the combustor wall. Predicted product con- 
centrations are nearly uniform across the combustor at this 
axial location. The experimental data shows considerable 
variability. Carbon dioxide concentrations vary more than 
the water. Product concentrations are overpredicted. The 
hybrid pdf simulation best predicts product concentration. 

Experimental data showed nearly complete reaction at an 
axial distance of 137.5 cm. The complete length of the 
combustor was not simulated. The actual combustor was about 
50% longer than simulated. In the pdf simulations, the much 
smaller distances in the radial direction cause the number 
of particles exchanged in the radial direction to vastly 
outnumber the number of particles exchanged in the axial di- 
rection. This hampers information transfer in the axial 
direction. Increasing the number of particles increases the 
number of particles exchanged axially, but is more computa- 
tionally expensive. The number of particles used in this 
■calculation was 250 particles per cell. Using low numbers of 
particles resulted in a gradual prediction of flame or 
combustor blow-off. Using fewer particles caused larger 
temperature oscillations in the iteration process. 

Particles from three different computational cells for 
an axial distance of 9.18 cm from the combustor step are 
shown in figure 5.78. The particles are plotted as a func- 
tion of mixture fraction and temperature, where mixture 
fraction is the sum of the mass fraction of fuel in all 
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Figure 5.78 Particle plots at an axial distance of 9.18 cm 
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species. Pure mixing with no combustion is represented by a 
line between mixture fraction of zero and temperature of 589 
K (air) to a mixture fraction of 1.0 at a temperature of 300 
K, which is pure fuel. The line would have some slight 
curvature due to variable C p . Combustion is represented by 
an elevation off this line. Complete combustion is repre- 
sented by a line from a mixture fraction of zero and temper- 
ature of 589 K (pure air) to a stoichiometric mixture frac- 
tion of 0.055 at 2500 K and hence to a mixture fraction of 

1.0 at a temperature of 300 K. 

The particle plots show that a range of mixtures is 
being simulated in various cells. The particle plot for the 
cell closest to the fuel inlet (r=1.35) shows various combi- 
nations of air and fuel with very little combustion. The 
maximum mixture fraction in this cell is 0.55. No pure fuel 
particles are shown. The plot for a radial location of 2.78 
cm shows very lean, but complete combustion. The particle 
plot for a radial distance of 5.62 cm shows slightly richer 
combustion and higher temperatures. 

Particle plots for an axial distance of 39.6 cm are 
shown in figure 5.79. These particle plots show much higher 
temperatures. As before, the level of combustion increases 
with distance from the combustor axis. The cell with a 
radial distance of 1.35 cm displays the largest number of 
unburned and partially burned particles. The cell at a 
radial distance of 5.62 cm displays only a few partially 
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Figure 5.79 Particle plots at an axial distance of 39.6 cm 
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burned particles. The reaction model implemented in the 
hybrid pdf model gives partially combusted particles for low 
temperatures, but most high temperature particles are almost 
fully burnt. 

The use of a finite-rate reaction scheme here has al- 
lowed significant overlap of fuel and oxidant concentra- 
tions. The assumed shape pdf, equilibrium chemistry model 
used by Smith and Smoot predicted little overlap of oxygen 
and fuel concentration profiles. The finite-rate combustion 
models used here underpredicted the carbon dioxide concen- 
trations in the recirculation zone, but the temperature 
predictions were much improved from their calculation. The 
hybrid pdf initially overpredicted fuel concentration, but 
predictions using the hybrid pdf were better than the base- 
line predictions at successive axial stations. Hybrid pdf 
predictions of oxygen, carbon dioxide and water were almost 
always superior to the baseline simulation. 
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5.3.2 SWIRLING DIFFUSION FLAME 

Attempts were made to simulate the combustor of Jones and 
Wilhelmi 121 . This is a newer set of experimental data with 
LDA measurements of axial, radial, and tangential velocity, 
temperature and species data at numerous axial stations 
including data very close to the injector face. This com- 
bustor flow had a 45 degree gas turbine type swirler, was 
100 mm in diameter and was 300 mm long. The fuel was gas- 
eous propane introduced by a central conical annular injec- 
tor. The axial velocity profile was completely positive by 
100 mm, but the large centerline defect was still apparent 
at 289 mm, near the combustor exit. All SIMPLE calculations 
showed a nearly uniform axial velocity at this position. 

This axial velocity defect may be caused by unsteady ef- 
fects. Also, there was some ambiguity in the initial axial 
velocity measurements. The first measurements showed two 
axial velocity peaks with a single radial velocity peak. 

This implied that the secondary peak was from the recircu- 
lating flow. This recirculating flow was stronger than any 
calculated. The second velocity peak may have something to 
do with a quarl at the top of the inlet swirler. The simu- 
lations did not model the quarl. At the second measurement 
station, the axial velocity exhibited strong reverse flow 
next to the fuel and air inlet flow, tapering off towards 
the centerline. All simulations attempted calculated the 
largest negative axial velocity along the centerline. The 
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experimental flow along the combustor centerline was quies 
cent. The high reverse flow next to the inlet flow and 
adjacent quiescent flow seems more characteristic of flow 
being entrained into a jet. 

It was decided to study a similar combustor with weaker 
swirl, as it is known that the k-e turbulence model poorly 
predicts areas of high swirl. A similar combustor is that 
of Jones and Tober 122 . This combustor has a thirty degree 
swirler, with no quarl. The fuel is propane which enters 
through a conical annular nozzle with a cone angle of 90 
degrees just inside of the air swirler. The internal diame- 
ter and length of the combustor are 196 mm and 700 mm, 
respectively. The inner and outer diameters of the swirler 
are 23.5 mm and 42 mm. The combustor was run at two air to 
fuel ratios. The case studied in this work had an air flow 
of 26.7 g/s and fuel flow of 1.37 g/s. 

A 70 by 55 grid rectilinear grid system was used. Grid 
stretching was used to concentrate grid points near the 
inlet fuel injector. The axial inlet velocity was calculat- 
ed by projecting back the velocity data from the x/D=.l data 
using the flow angle from the peak axial and radial veloci- 
ties, and then scaling the profile to give the correct mass 
flow. The radial velocity inlet data was found by doing nu- 
merous runs until a good match was obtained for the experi- 
mental data at the first couple of measurement stations. 

The hybrid pdf simulation used a stronger inlet radial velocity. 
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The flowfield is characterized by two hot recirculating 
flows. One recirculating flow is a toroidal vortex adjacent 
to the combustor walls. This recirculation zone extends to 
about x/D=0 . 5 . The second recirculating flow occurs along 
the combustor axis and is referred to as a Central Recircu- 
lating Zone or CRZ. This recirculation zone is characteris- 
tically caused by high inlet swirl. Based on a contour plot 
in Jones and Tober, the CRZ appears to extend to x/D=1.6, 
and takes up half the diameter of the combustor at x/D=0.8. 
The combustor was simulated out to x/D=2.5. The two recir- 
culation zones are separated by higher velocity flow which 
is at an angle of about 30 degrees to the combustor axis. 
Initially, this jet-like flow is largely composed of air and 
fuel from the combustor inlet. This flow is typically fuel 
rich at smaller radii and oxygen rich at larger radii. 

Predictions and measurements of velocity and tempera- 
ture at the x/D= . 1 location are shown in figure 5.80. The 
calculated axial velocity peak for the flow separating the 
two recirculation zones is well predicted. This was one of 
the criteria for selecting inlet conditions. The width of 
the peak is overpredicted. The pdf simulation better pre- 
dicts the lower portion of the velocity peak. The experi- 
mental data exhibits nearly constant, low velocity flow 
across the CRZ. Both calculations show the magnitude of the 
reverse-flow velocity increasing towards the combustor axis. 
The predicted maximum reverse flow velocity in the CRZ is 
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Figure 5.80 Velocity and Temperature at x/D = 0.1 
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much higher than measured, which was characteristic of trial 
calculations of Jones and Wilhelmi's combustor. The outer 
recirculation zone has rather low velocity which is well 
predicted by the simulations. Near the outer combustor 
wall, calculations exhibited slightly larger negative axial 
velocity than measured. 

Radial velocity predictions for the outer recirculation 
zone are very good for both simulations. Both simulations 
underpredict the radial velocity peak and overpredict the 
width of the higher velocity jet-like flow by almost a 
factor of two. Measurements of radial velocity across the 
CRZ are nearly zero. The simulations predict slightly 
higher radial velocity flow. The radial velocity profiles 
near the combustor axis are parabolic, and scale on the 
cross-sectional size of the CRZ. The hybrid pdf radial 
velocity profile for the jet-like inlet flow is slightly 
shifted outward due to a larger inlet radial velocity pro- 
file. In this and other simulations of this flowfield, the 
pdf simulations exhibited lower radial displacement of the 
inlet flow with axial distance. 

The predicted tangential velocity peaks are under-pre- 
dicted and are highly diffusive. The experimental measure- 
ments show very low tangential velocity in the CRZ. The 
simulations show nearly constant angular rotation across the 
CRZ and up to the location of peak tangential velocity. 
Tangential velocity in the outer recirculation zone is well 
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predicted by both simulations. The simulations do not 
predict an experimentally measured depression in tangential 

velocity at a radius of 7 cm. 

The numerical simulations correctly predict high temper- 
atures in both of the recirculating zones. However, these 
temperatures are over-predicted by hundreds of degrees. The 
high temperature predictions are the result of assuming 100 
% combustion efficiency. Temperature predictions would im- 
prove if a lower heat release were used. The pdf module 
used here uses thermodynamic data of Gordon and McBride and 
assumes 100 % efficiency. This model is more difficult to 
manipulate to predict lower heat of reaction. The experi- 
mental data shows depressed temperatures near the outer 
wall. The combustor wall was treated as adiabatic in the 
both of the simulations as the pdf module does not have a 
wall heat transfer model in it. Thus, the predictions do not 
exhibit a drop in temperature near the wall. The large drop 
in temperature between the two recirculation zones is pre- 
dicted, but the minimum temperature predicted is hundreds of 
degrees too high. Both simulations predict the same minimum 
temperature . 

Fuel, oxygen, carbon dioxide, and mixture fraction 
concentrations at x/D=.2 are shown in figure 5.81. The 
mixture fraction is a combination of the mass fractions of 
hydrogen and carbon in all of the species. Other concentra- 
tions are in mole fractions. Mixture fraction, Carbon 








207 


dioxide predictions and measurements agree in the outer 
recirculation zone. The traditional type simulation over- 
predicts carbon dioxide concentration at other locations. 

The conventional simulation overpredicts the minimum carbon 
dioxide concentration by a factor of almost three. The 
hybrid pdf simulation does a better job of predicting carbon 
dioxide concentration. Both simulations show carbon dioxide 
initially increasing with radial distance from centerline. 

The experimental concentrations are constant in this area. 

The mixture fraction and carbon dioxide profiles match in 
the outer recirculation zone as the fuel is fully burned in 
this region. The drop in mixture fraction at a radius of 4 
cm isn't shown in the simulations. This is due to the 
overprediction of carbon dioxide at this radial location. 

The peak experimental mixture fraction occurs at a larger 
radius than the peak for carbon dioxide. Thus, there must 
be significant fuel at a radius of 3 cm. The simulations 
•predict significant fuel at this location. Unfortunately, 
the predicted fuel concentration is not large enough to 
predict the mixture fraction peak at the correct radial 
location. Instead, the simulations predict mixture fraction 
peaks at a lower radius where the fuel and carbon dioxide 
concentrations overlap. The conventional simulation pre- 
dicts the magnitude of peak mixture fraction. Unfortunate- 
ly, this is caused by an over prediction of carbon dioxide. 
Curiously, the conventional simulation predicts peak fuel 
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concentration in the same region as peak carbon dioxide. 

The oxygen concentration predictions show the same trends as 
the experimental data. The predicted concentration of 
oxygen in the outer recirculation zone is excellent, as was 
other concentrations. At a radius of 4 cm, a peak is shown 
due to oxygen from the inlet flow. At 2.5 cm, the oxygen 
concentration is at a minimum due to large concentrations of 
carbon dioxide and fuel. The conventional simulation under- 
predicts oxygen, especially across the inlet flow. The pdf 
simulation best predicts the oxygen concentration. The 
radial location of peak oxygen concentration is well pre- 
dicted. The hybrid pdf predicted peak lies midway between 
the experimentally measured peak and conventionally predict 
ed peak. The low oxygen concentration inside the CRZ is 
best predicted by the pdf simulation, but the width of this 
part of the profile is predicted by the conventional simula- 
tion. 

Temperatures corresponding to the above concentrations, 
are shown in figure 5.82. Temperature predictions in the 
recirculation zones are high. The temperature predictions 
across the outer recirculation zone are nearly constant, as 
were the species concentrations. The minimum temperature 
predicted by the pdf simulation is midway between the exper- 
imental minimum and conventionally predicted minimum. At 
the previous axial location, the simulations predicted the 
same minimum temperatures. The radial location of minimum 
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Figure 5.82 Temperature at x/D = .2 
temperature is the same as that for minimum carbon dioxide 
and maximum oxygen concentration. The hybrid pdf simulation 
best predicts temperature in the CRZ. 

Velocity and temperature at x/D=.3 are shown in figure 
5.83. Both simulations predicted nearly identical axial 
velocity profiles. Predicted positive axial velocity pro- 
files are much wider than the experimental profile. Peak 
axial velocity is underpredicted, although slightly beyond 
the peak, the predictions are excellent. The size of the 
CRZ is underpredicted. The maximum reverse velocity is 
over-predicted in the CRZ, although it is much improved from 
the x/D= . 1 location. Again, the experimental measurements 
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show nearly constant, low axial velocity across the CRZ. 

The experimental radial velocity profile shows a narrow high 
velocity peak. Predicted peak profiles are blunter and 
wider. The general shape of the predicted profiles is very 
similar, only differing in magnitude. The hybrid pdf simu- 
lation predict a velocity profile about 30 % lower than the 
conventional simulation. 

The experimental tangential velocity profile shows two 
peaks. The first peak is due to the swirling inlet swirling 
flow. The second velocity peak is adjacent to the outside 
combustor wall. Tangential velocities in this region are 4- 
6 m/s. At x/D=0.1, tangential velocities were less than 5 
m/s near the outer wall. The increase in tangential velocity 
at larger radii is due to spreading of the peak tangential 
velocity which has decreased by about a factor of two from 
the x/D=0 . 1 location. The numerical simulations only pre- 
dict the first velocity peak. Examination of the axial and 
radial velocity profiles does not give an apparent reason 
for the second outer tangential velocity peak or a tangen- 
tial velocity defect in the profile. 

The shape of the predicted temperature profiles is 
qualitatively correct. Predicted temperatures are too high. 
The conventional simulation predicts the minimum temperature 
at about the correct radial location. The hybrid pdf pro 
file shows the radial location of minimum temperature about 
one cm less than was measured. Both simulations signifi- 
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cantly underpredict the radial position of peak temperature. 
The pdf simulation predicts lower temperatures than the 
conventional simulation. Experimental temperature gradients 
are much larger than predicted. 

Velocity and temperature results at x/D=.4 are shown in 
figure 5.84. The agreement between the simulations is 
better than the agreement between the predictions and mea- 
surements . The simulations predict the same magnitude of 
peak axial velocity, but the hybrid pdf simulation predicts 
a smaller radial displacement of the velocity profiles. 
Predicted profile peaks are again wider and blunter than the 
experimental peaks. Predicted radial velocities are almost 
identical out to a radius of 4 cm. Beyond 4 cm, the hybrid 
pdf simulation predicts lower radial velocities. The veloci- 
ty predictions within the CRZ are closer to the experimental 
measurements than at previous stations. The width of the 
CRZ is again underpredicted. The maximum reverse flow 
velocity in the CRZ is decreasing. The experimental velocity 
gradient between the reverse flow in the CRZ and adjoining 
flow is greater than predicted. Predicted velocity profiles 
show diffusive, gradual changes. The comparison between 
predicted and measured tangential velocities is very good 
out to a radius of 3 cm. Beyond 4 cm the hybrid pdf simula- 
tion predicts higher tangential velocity at a lower radii, 
which agrees with the data. This is due to the reduced 
radial expansion predicted in the hybrid pdf simulation. 
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Again, experimental data shows a second tangential velocity 
peak, which isn't predicted by the calculations. Temper- 
ature predictions are high. Temperatures continue to become 
more uniform with axial distance. The pdf simulation pre- 
dicts lower temperatures than the conventional simulation, 
but the conventional simulation better predicts the radial 
location of peak temperature. 

Concentrations for x/D=.4 are shown in figure 5.85. 

The simulations predict relatively constant carbon dioxide 
concentration across the combustor. The experimental mea- 
surements predict a 30 % drop in concentration at a radius 
of 7 cm. The conventional simulation predicts peak carbon 
dioxide at 5 cm. Carbon dioxide is underpredicted at small- 
er diameters and overpredicted at larger diameters. The 
experimental mixture fraction shows only a slight dip at 7 
cm, unlike the carbon dioxide profile. Also, the mixture 
fraction shows a significant peak at 5.5 cm while the peak 
carbon dioxide is at 3 cm. This implies that there is 
significant fuel or partially burned fuel between these 
points. The fuel mole fraction should be on the order of 
0.01. The conventional simulation predicts almost no fuel 
at this axial station. The pdf simulation shows a maximum 
mole fraction of 0.0025 around 5 cm. 

The conventional simulation does a good job of predict- 
ing oxygen concentration out to 6 cm. The peak experimen- 
tal oxygen concentration at 7 cm isn't predicted in the 
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simulations . The pdf simulation predicts the highest oxygen 
concentration at this station. Neither simulation is supe- 
rior in the prediction of oxygen. Both simulations overpre- 
dict combustion. This results in more expansion at smaller 
axial and radial displacement from the combustor inlet. 

Velocities at x/D=.6 are shown in figure 5.86. Nega- 
tive axial velocity is measured over 60 % of the combustor 
inner diameter at this station. The simulations underpredict 
the size of this flow. Thus, the experimental velocity 
profiles show higher positive axial flow at greater radii. 
The back-flow velocity along the combustor axis is well 
predicted. The hybrid pdf simulation shows lower reverse 
velocity near the combustor axis, which may be due to in- 
creased combustion at smaller radii. The conventional simu- 
lation best predicts the axial velocity profile at this 
location. Experimental radial velocities are all positive. 
Simulations predict negative radial velocity for at least 
half of the inner diameter. The pdf simulation shows higher 
radial velocity at larger radii, unlike previous stations. 
Both simulations predict very similar tangential velocity 
profiles at this station. Predicted tangential velocities 
are nearly constant from 3 cm to 9 cm. The experimental 
data shows a tangential velocity peak at 5 cm, which is 
significantly underpredicted in the simulations. 

Species concentration and temperature at x/d=.6 are 
shown in figure 5.87. The experimental carbon dioxide 
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concentration is constant out to a radius of 8 cm. The 
mixture fraction shows a slight peak at 8 cm, which implies 
that there is still unreacted or partially reacted fuel at 
this location. Insufficient fuel is predicted at this axial 
station. The shape of the predicted mixture fraction pro- 
file is the same as the carbon dioxide profile. Measure- 
ments show a temperature peak and a drop in oxygen concen- 
tration at 8 cm, which suggests ongoing combustion at that 
location. The simulations show increased oxygen concentra- 
tion past 6 cm. 

Velocity and temperature at x/D=0.8 are shown in figure 
5.88. The experimental data show the CRZ is about the same 
diameter as the previous station, but the reverse flow 
velocity has slightly increased. The simulations show a 
smaller diameter CRZ and reduced reverse flow velocity from 
the previous station. Predicted radial velocities are all 
negative at this station. The experimental data for radial 
'velocity is rather sparse, but suggest a different flow 
pattern. The discrepancy in prediction of peak tangential 
velocity grows to 50 %. The second tangential velocity peak 
isn't seen at this location. Measured tangential velocities 
at this station are greater than at the previous axial 
location at all radii. This implies there is some error in 
tangential velocity. Predicted temperatures across the 
combustor are nearly constant. The experimental data shows 
a slight temperature peak at 9 cm. The pdf simulation shows 
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Figure 5.88 Velocity and Temperature at x/D=0.8 
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lowest temperature, closer to experimental data. 

Mixture fraction isn't perfectly conserved in these 
calculations. If everything were perfectly conserved in the 
calculations, the predicted temperatures would be closer at 
this location. Summing the conventional mixture fraction 
flow across the combustor shows a slight increase in total 
mixture fraction at this location. Summing the hybrid pdf 
predictions of mixture fraction flow shows a slight de- 
crease. Thus, some of the improvement seen between the 
calculations should be attributed to error, not better 
modeling. 

Some cells on which to perform particle plots were 
chosen for examination. Particle plots of mixture fraction 
versus temperature for x/D of 0.0034 are shown in figure 
5.89. The particle plot for r/D of 0.054 displays a wide 
variety in mixture fraction and temperature. This cell is 
within the CRZ . Most particles that appear in this cell are 
fully reacted. The particle plot for r/D of 0.74 displays 
three distinct particle groupings. One grouping is air with 
a slight amount of unburned fuel. The second group is a 
single particle at 1914 K and a mixture fraction of 0.14. 

The third particle grouping is almost pure unburned fuel. 
This particle plot is of a cell between the two recircula- 
tion cells. The particle plot for r/D of 0.198 is actually 
of 250 particles clustered around 1940 K and a mixture frac- 
tion of 0.044. The fluid from this cell is within the outer 
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Figure 5.89 Particle plots at x/D=0.0034 
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recirculation zone, next to the combustor wall. The fluid 
is almost completely mixed. The fluid velocity in this cell 
is very low. Information transfer in this zone of the 
combustor is very slow in the hybrid pdf calculation. 

A second set of particle plots for x/D of 0.136 is 
shown in figure 5.90. All of these plots show considerable 
variability. No particularly fuel rich particles are shown. 
The plot for r/D=0.054 is within the CRZ. Most particles 
are fully burned except for a smaller number centered around 
900 K. The plot for r/D=0.074 is closer to the fuel stream. 
The number of partially burned particles centered near 900 K 
is larger than for the r/D=0.054 plot. There are also a few 
more particles with higher mixture fraction. The plot for 
r/D= . 198 has the greatest number of particles with moderate 
temperature and mixture fraction. These particle plots are 
used to show variability within the fluid cells. Particles 
can be plotted over one another. Some sort of statistical 
analysis would probably help in interpreting this data. 
However, much of the easier statistical analysis assumes 
normal or Gaussian type distributions. The distributions 
are obviously bimodal at low temperature in and near the low 
temperature fuel-air stream. 

Neither simulation captures the experimentally mea- 
sured low fluid velocity of the CRZ near the combustor 
inlet. The pdf simulation predicts slower combustion and 
generally less radial displacement. The inlet radial 
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velocities had to be increased so that the hybrid pdf re- 
sults would compare to the conventional simulation. While 
the lower rate of combustion predicted by the pdf simulation 
predicted temperatures closer to the experimental data, the 
larger radial displacement of fuel and its products by the 
conventional simulation produces profiles that were closer 
to the data in regards to shape. Both simulations use 
hybrid differencing for calculation of coefficients. As the 
flow is about 30 degrees to the grid (initially 45 degrees 
for fuel) , hybrid differencing can be expected to give 
larger errors due to numerical diffusion, compared to higher 
order schemes . The pdf scheme should be tested with higher 
order schemes. Unfortunately these schemes will to have all 
positive coefficients, if the same differencing scheme is 
used to calculate the pdf particle number exchanges as 
transporting a negative number of particles does not have a 
physical meaning. A quick fix for this may be to use lower 
order differencing for the transport of particles. 
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5.3.3 PREMIXED COMBUSTION CASE 

To explore the further applicability of the hybrid pdf 
model, a lean premixed combustor case is considered. In a 
premixed combustor, fuel and oxidant are mixed upstream of 
the actual combustor. This results in uniform air to fuel 
ratios. In normal combustors, there is a wide variation in 
air to fuel mixture. Some areas are lean, some are rich, 
and some are near stoichiometric. Burning near stoichiomet- 
ric results in high temperature which markedly increases 
nitrous oxide emissions. Lean premixed combustors do not 
have this problem. Thus, lean premixed combustion is one 
way in which to lower nitrous oxide emissions. Nitrous 
oxides are one component of smog. Currently, some European 
countries are proposing taxes based on nitrous oxide emis- 
sions of the aircraft flying through their airspace. Since 
the combustor inlet mixture is premixed and is combustible, 
the location where combustion occurs must be fixed or an- 
chored. Combustion can flash back into the inlet duct or 
can be blown downstream. This is deleterious to performance 
and safety, particularly in aircraft. Combustion is an- 
chored by forming hot recirculation zones. One way to form 
a recirculation zone is to put an obstruction, or flame 
holder, in the combustor. Some examples of flame holders 
are disks, disks with holes, v-gutters, and conical bluff 
bodies. In the combustor studied here, a large increase in 
flow area is used to form a recirculation zone. This con- 


MT! 



227 


figuration is also known as a step or dump combustor. The 
basic flow configuration is the familiar flow over a step. 

The combustor studied here is that of Gould, Stevenson, and 
Thompson 123 at Purdue University. Premixed fuel and air 
goes through a converging nozzle, with an exit diameter of 
76.2 mm, into the combustor chamber of 152.4 mm outer diame- 
ter. The inlet nozzle was contoured to provide a nearly 
uniform velocity flow across the combustion chamber inlet 
plane. There is a slight velocity defect due to the bound- 
a.ry layer. The inlet mixture is gaseous propane and air at 
an equivalence ratio of 0.5. Velocities were measured using 
Laser Doppler Velocimetry. Temperature was measured using 
high speed thermometry. Various turbulent correlations were 
also measured. A corrective lens was used to take optical 
measurements. This lens degraded due to combustion, limit- 
ing the number of axial stations where data was taken. 

Gould performed computational combustion simulations as 
part of his dissertation study. A pressure based solver 
with the eddy-breakup combustion model of Magnussen and 
Hertajer 124 was used. In this combustion model, combustion 
is proportional to turbulent eddy lifetime, and minimum 
species concentration of fuel, oxygen or combustion product. 
Combustion is initiated in this model by temporarily intro- 
ducing product species at specified cells in the simulation. 
Recently Magnussen 125 proposed an improved combustion model 
called the Eddy Dissipation Concept combustion model. This 
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model was used in a calculation of Brum and Samuelson's 
combustor 126 , which was then compared to the predictions of 
Nikjooy [28] . 

An inlet axial velocity of 22 m/s was specified in the 
combustor calculations, except for the boundary layer. Zero 
radial velocity was specified. There was no swirl compo- 
nent. The measured inlet turbulent kinetic energy was very 
low, 0 . 0006*U 2 , except for the boundary layer, which was 
measured at 10.0 m 2 /s 2 . Combustion simulations did not 
satisfactorily converge using such low inlet turbulence. 
Values of 0.0054*U 2 were used, except for the boundary 
layer. This allowed the simulations to converge much fur- 
ther. High inlet temperature was initially specified to 
initiate combustion simulations. The simulations were then 
restarted using the correct inlet temperature of 300 K. The 
pdf simulation started from a converged SIMPLE calculation 
restart file. As the current pdf model does not allow heat 
transfer to the wall, adiabatic wall treatment was also used 
in the conventional simulation. The heat of reaction term 
used in the conventional simulation was adjusted to give the 
same final temperature as the thermodynamic model used in 
the pdf module for a single cell calculation. 

Experimental values of velocity, temperature, and 
various double and triple correlations of these quantities 
were measured. The experimental data taken near the combus- 
tor axis shows radial velocities around -2.0 m/s. This is 
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physically unrealistic. Evidently, the thermometry appara 
tus affected radial velocity measurements. Radial veloci- 
ties measured without the thermometry apparatus for the 
isothermal case were reasonable. Predicted radial velocity 
profiles do not agree with the measured profiles for the 
burning case. There was significant anisotrophy in the 
turbulent Reynold's stresses. Axial velocity and tempera- 
ture at a location of one step height past the combustor 
inlet are shown in figure 5.91. The axial velocity is well 
predicted by the two different combustion model simulations. 
The axial velocity in the main flow is slightly over-pre- 
dicted. The experimental data shows a slightly sharper 
velocity gradient between the cold core flow and the hot 

I 

recirculation zone. Temperature measurements in the recir- 
culation zone are hundreds of degrees lower than predicted 
temperatures. Low temperatures are measured near the com- 
bustor wall. The predictions assumed no heat transfer to 
the walls. Both the conventional simulation and hybrid pdf 
simulation predict complete combustion for much of the 
recirculation zone. As velocity is lowest next to the 
combustor wall, this portion of the flow has the longest 
time for combustion to occur. With the high temperatures 
measured at other axial stations, much of the recirculation 
zone should be fully reacted. This researcher feels that 
the peak experimental temperature at this axial location 
should closely match the peak temperatures found at other 




OOO 





231 


axial stations. Allowing some heat transfer in the calcula 
tions would improve the temperature predictions along the 
wall, but shouldn't produce a 400 degree drop in temperature 
across the whole recirculation zone. Experimental velocity 
data was only taken out to 80 % of the combustor diameter. 
Predictions of axial velocity in the recirculation zone 
appear quite good, despite the large temperature discrepancy 
in the recirculation zone. This implies that temperature 
measurement in the recirculation zone at this station is in 
error or, velocity is not as strongly coupled to density as 
it is in the core flow. 

The results at an axial location of 3 step heights are 
shown in figure 5 . 92 . The axial velocity predictions and 
measurements are very close. The pdf predictions show a 
slightly higher axial velocity outside of the cold central 
core. The hybrid pdf simulation shows a slightly faster 
combustion progress into the core flow. High temperature is 
predicted at slightly lower radii in the pdf simulation. 

The pdf simulation predicts a wider blunter temperature peak 
than the conventional simulation. The conventional simula- 
tion predicts slightly higher temperature than measured in 
the hot recirculation zone. 

At 5 steps heights past the inlet, the axial velocity 
predictions, shown in figure 5.93, are very similar. The 
conventional simulation comes closest to the measurements in 
the outer recirculation zone, but axial velocity is 
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bracketed by the predictions in the inner recirculation 
zone. The small negative axial velocities shown near the 
wall in the conventional simulation mean that the recircula- 
tion bubble is almost closed for that simulation. The pdf 
simulation shows slightly positive axial velocity near the 
wall. The axial velocity near the combustor axis is under- 
predicted by both the simulations. Previously velocity was 
slightly overpredicted. In this same area, predictions and 
measured data for temperature are very close. Beyond the 
cold core flow, the hybrid pdf simulation noticeably pre- 
dicts higher velocity and temperature than the both the 
conventional simulation and experimental measurements. Up 
to a radius of 4.5 cm, the conventional simulation predicts 
temperature almost perfectly. At the previous station, 
temperature was slightly overpredicted. Beyond a radius of 
4.5 cm, temperature measurements are bracketed on the low 
side by the conventional simulation and on the high side by 
'the hybrid pdf simulation, except near the wall. Qualita- 
tively, the shape of the pdf simulation temperature profile 
is closest to the measured temperature profile at this axial 
location. 

At the last measurement station of 12 step heights past 
the inlet, both numerical simulations deteriorate, as shown 
in figure 5.94. The difference between the simulations and 
experimental data is largest at this station. In most 
experimental combustors the last measurement station is 
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usually at the completion of combustion and mixing, and exit 
profiles are fairly uniform. In this experiment, a large 
fraction of the cold inner flow is still unreacted. Thus, 
there are significant gradients in combustor profiles. Both 
predictions show a higher axial velocity along the center- 
line than measured, while temperature is underpredicted. 

This is somewhat of an anomaly. Beyond a radius of 2 cm, 
the conventional simulation underpredicts axial velocity, 
while the pdf simulation continues to overpredict it . The 
experimental measurements show an almost linear increase in 
temperature off the combustor axis. At a radius of two 
centimeters, both temperature predictions are off by 300 
degrees K. The hybrid pdf simulation best predicts the high 
temperature zone at this temperature. Gould noted similar 
deficiencies in his predictions and theorized that it was 
the fault of turbulence model. At previous stations, the 
predictions in and near the low turbulence level cold core 
flow were excellent. At this last station, the cold flow is 
hotter and more turbulent. The results at this station 
could be more closely predicted by increasing inlet turbu- 
lence. Unfortunately, this would work at the detriment of 
previous predictions. With the extremely low inlet turbu- 
lence, turbulence modeling may be more in the realm of 
turbulent transition. The turbulence model used in these 
simulations assumes high levels of turbulence, that is, it 
is a high-Reynolds number turbulence model. It would be 
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interesting to try a low Reynold's number turbulence model 
in another simulation. 

The hybrid pdf temperature simulation exhibited sharp- 
er temperature gradients between the cold inner flow and hot 
outer flow, and higher temperature, particularly at the last 
station. At the last station, the pdf simulation best pre- 
dicts temperature. 
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5.3.4 SWIRLING HYDROGEN DIFFUSION FLAME 

Recently, a group of researchers at the University of 
Dayton, Research Institute, have taken benchmark data for 
model validation. This includes non-swirling and swirling 
air jet flows 127 , methane jet flames 128 , and hydrogen jet 
diffusion flames 129 . The combustion data for the burning 
hydrogen is the most complete. Three components of velocity 
were measured along with many turbulence quantities using 
LDV equipment. Temperature was measured using Coherent 
Anti-Stokes Raman Spectroscopy (CARS) . Typically 500 tem- 
perature measurements were taken at each location. These 
measurements were processed to get rms temperature fluctua- 
tions and temperature pdfs. The CARS method is still being 
developed and the data is ideal for validating the composit- 
ion pdf method employed here. 

The experimental system consists of a central fuel tube 
of 9.45 inner diameter, a concentric annular air tube of 
26.92 inner diameter, and external non-swirling air supplied 
to a vertical 150 by 150 mm semi -rectangular test section 
which is 486 mm long. Various combinations of velocity and 
swirl were tested. The case studied here is case no. 3 
which had a fuel jet bulk velocity of 100 m/s, annulus air 
bulk velocity of 20 m/s with thirty degree swirl, and exter- 
nal air bulk velocity of 4 m/s. Radial data was taken at 
1.5, 10, 25, 50, 150, 250 mm from the inlet. Additional 
combustor centerline data were also taken. Data at the 
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first measurement plane were used for inlet data in the 
simulations. The mass fraction of species were inferred 
from the geometry and temperature measurements. The high 
temperature combustion interface between the fuel and air 
jets provided an ignition source. A one step chemistry 
submodel was employed. This mechanism was inferred from 
several larger reaction mechanisms for heavier fuels which 
included a hydrogen-oxygen to water reaction step. The 
activation energy used was 40 Kcal/mole. The exponent used 
for the fuel concentration was 0.25, and the oxygen exponent 
was 1.5. The frequency number had to be increased to 2.7 x 
10 13 to simulate continuing combustion in the calculations. 

Baseline calculations were done with a coarse grid of 
40 by 30 and a fine grid of 70 by 60. These simulations 
predicted nearly identical axial velocity and temperature 
profiles. The pdf calculations were performed using the 
same grids. An adiabatic wall boundary condition was im- 
posed at a radius of 75 mm. 

Temperature, turbulent kinetic energy, axial and tan- 
gential velocity at an axial distance of 25 mm from the 
inlet plane are shown in figure 5.95. Temperature at this 
plane is best predicted by the baseline coarse grid simula- 
tion. The coarse grid temperature profiles are largely the 
same as the fine grid profiles except the coarser grids 
truncate the temperature peaks. The peak temperature is 
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slightly overpredicted by the baseline fine grid calcula 
tion. The pdf simulations underpredict temperature. Ef 
forts to increase the coarse grid peak temperature by in- 
creasing reaction rate and maximum inlet temperature at this 
axial location were largely unsuccessful. 

There are three sets of experimental measurements. 
Measurements of velocity and turbulence quantities are 
conditional. That is, the fuel jet, annular jet, and exter- 
nal air were seeded in successive runs. Calculated turbu- 
lence kinetic energy and velocities are not conditional. An 
average, unconditional velocity is calculated. To calculate 
conditional velocities a joint velocity-scalar pdf would 
have to be developed. 

Predictions of turbulent kinetic energy are very good. 
The baseline fine grid calculation best predicts the peak 
turbulent kinetic energy. The fine grid pdf simulation 
predicts slightly higher peak turbulent kinetic energy. The 
two coarse grid solutions predict the highest peaks. The 
pdf simulations best predict centerline and near centerline 
turbulence kinetic energy. At larger radii, the fine grid 
calculations slightly overpredict the turbulence kinetic 

energy. 

Axial velocity predictions tightly bracket the exper- 
imental measurements out to a radius of 5 mm. The baseline 
simulations predict higher axial velocity due to the higher 
predicted temperatures. The fine grid pdf simulation best 
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predicts the axial velocity deficit between the fuel jet and 
annular air flow. At larger radii, predicted axial veloc- 
ities converge. The pdf simulations predict slightly 
higher tangential velocity out to a radius of 12 min due to 
lower temperature predictions and reduced expansion of the 
flow. The tangential velocity predictions are quite good. 
The coarse grid calculations predict higher tangential 
velocity beyond the tangential velocity peak more in agree- 
ment with the data. Overall, the baseline coarse grid 
calculation best matches tangential velocity data. 

Some plots of temperature versus mixture fraction 
(total fuel component) at an axial distance of 26 mm are 
shown in figure 5.96 for the coarse grid solution. There 
are 250 particles in each cell plot. Each computational 
grid cell in this region is of the order of one mm in the 
radial direction and 3 mm in the axial direction. The cells 
are in the mixing layer region where the combustion rate is 
highest. The pdf simulation predicts a wide variation of 
temperature and mixture fraction in the mixing layer. Each 
particle generally falls close to a state function composed 
of a straight line connecting mixture fraction of zero, 
temperature 300 degrees (pure air) with mixture fraction 
0.029, temperature 2550 (stoichiometric or complete burning 
of both reactants) and thence a curved line ending at mix- 
ture fraction of 1.0, temperature of 300.0 (pure fuel). 
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This composite line represents complete burning of fuel for 
lean mixtures and complete consumption of oxidant for rich 
mixtures, effects of variable C p , and a fixed ratio of 
nitrogen to oxygen. Temperature is highly nonlinear with 
mixture fraction. Temperature increases most rapidly near 
stoichiometric. This thermodynamic behavior causes very- 
large temperature gradients in the mixing layer . 

The cell plots show that the total fuel component 
decreases with increasing radius. Only a few particles are 
richer than stoichiometric at a radius of 12.5 mm. There 
are at least a dozen partially reacted particles. These 
partially reacted particles have temperatures less than 850 
degrees. Most particles above 850 degrees are either fully 
reacted or very close to being fully reacted. As the major- 
ity of particles in these plots are fully reacted, this 
explains why efforts to increase reaction rate didn't in- 
crease peak temperature predictions. The mixing of fuel and 
oxidant species would have to be increased to increase 
temperature predictions in the pdf simulation. 

Initial coarse grid pdf simulations were done using 100 
and 250 particles per cell. Predictions of velocity and 
temperature were largely unaffected by varying the particle 
numbers. What did seem to affect initial calculations was 
the method of treating residual or partial particle trans- 
fers. The results shown here were obtained by keeping track 
of residual particle transfers and actually exchanging 
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particles when the residual was greater or equal to one, 
rather than randomly adding in particles as was done previ- 
ously. This improved temperature predictions at larger 
radii . 

Four cell plots of temperature vs mixture fraction for 
the fine grid pdf simulation are shown in figure 5.97. There 
are only 100 particles in these cells. Most of the parti- 
cles in these cells plots are fully reacted. 

Five coarse grid pdfs of temperature at an axial dis- 
tance of 26 mm are shown compared to experimental measure- 
ments at 25 mm in figure 5.98. The simulated pdf is taken 
from the final iteration. The temperatures passed to the 
velocity solver are averaged over hundreds of iterations. 
Thus , calculated pdfs should be also averaged over hundreds 
of iterations. This would require more memory than is 
currently needed for the complete calculation, unless a 
weighted running average is employed. Temperature measure- 
ments or particle temperature predictions are partitioned 
into bins which are multiples of one hundred degrees. The 
experimental measurements don't show any temperatures beyond 
2500 degrees while the simulations show a few particles just 
beyond 2500 degrees. The simulations assume 100 percent 
combustion efficiency. The experimental pdfs seem to be 
composed of one or two gaussian pdfs. At the largest radial 
location the predicted and measured temperature pdf are very 
similar, a single peak near 300 degrees. The pdfs taken at 
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the next smaller radial distance still show a pdf peak near 
300 degrees, but differ largely in the number of particles 
at higher temperatures. The simulation correctly predicts a 
wide variation in temperature at this location. At the 9.5- 
9.6 mm radial cell location, the experimental peak pdf 
occurs at higher temperature, while the simulation shows a 
pdf peak at much lower temperature. The average temperature 
plot showed an experimental peak temperature at 9.6 mm, 
while the coarse grid pdf simulation predicted a peak tem- 
perature at 8.3 mm. This agrees with the pdf plots shown 
here. The pdf coarse grid simulation shows the largest 
number of high temperature particles at 8.3 mm. The pdf is 
fairly uniform at this radius and this is one of the reasons 
why a fine grid pdf solution was done. 

A comparison of fine grid pdfs versus the same experi- 
mental pdfs is shown in figure 5.99. The fine grid pdf 
simulation predicts a peak temperature at 9 mm. This cell 
’isn't shown here, but the plot for 8.4 mm shows more parti- 
cles at higher temperature than does the coarse grid pdf at 
8 . 3 mm. 

Profiles of temperature, turbulent kinetic energy, 
axial and tangential velocity at 50 mm are shown in figure 
5.100. Peak temperature is best predicted in the coarse 
grid pdf simulation. At this station, the fine grid pdf 
simulation underpredicts peak temperature. The baseline 
simulations overpredict peak temperature by over 300 
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Figure 5.100 Temperature, turbulent kinetic energy, axial 
and tangential velocity at x = 50 mm. 
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simulations overpredict peak temperature by over 300 de 
grees. All numerical simulations predict very sharp temper- 
ature peaks while temperature measurements show more rounded 
peaks. The pdf simulations best predict turbulent kinetic 
energy. The baseline simulations predict larger radial 
displacement of the turbulent kinetic energy peak, overpre- 
dict the peak, and underpredict turbulence kinetic energy at 
and near the combustor axis. Axial velocity is best pre- 
dicted by the pdf simulations. Baseline simulations over- 
predict centerline axial velocity by about 20 %. The axial 
velocity predictions between the simulations converge with 
increasing radius. The baseline coarse grid simulation best 
predicts the tangential velocity profile. The coarse grid 
simulations are almost identical beyond a radius of 13 mm. 
The fine grid tangential velocity predictions are almost 
identical beyond a radius of 12 mm. The tangential velocity 
predictions are affected by predicted heat release (radial 
flow expansion) and also show grid dependence. 

Comparisons at an axial distance of 75 mm from the 
inlet are shown in figure 5.101. The pdf simulations under- 
predict the peak temperature by about 100 degrees, while the 
baseline simulations overpredict temperature by 300-400 
degrees. Predicted temperature peaks are quite sharp, while 
the temperature measurements are starting to show a plateau 
like area of high temperature just inside the peak tempera- 
The fine grid pdf simulation shows the bluntest 
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temperature peak. The pdf simulations excellently predict 
turbulent kinetic energy and axial velocity, while the 
baseline simulations overpredict both quantities out to a 
radius of 15 mm. From a radius of 9 to 15 mm, the predic- 
tions of axial velocity are surprisingly close to each 
other, despite the large difference in temperature predic- 
tions at this axial station. The largest difference in 
axial velocity predictions occurs along the flame center- 
line. The fine grid simulations predict the lowest tangen- 
tial velocity which is also closest to the experimental 
measurements. The coarse grid pdf simulation predicts 
highest tangential velocity, in part due to lower predicted 
heat release. At the previous axial station, the coarse 
grid simulations best predicted tangential velocity beyond 
the peak. Under a radius of 8 mm at this station, the shape 
of the tangential velocity appears to be poorly predicted. 
Predicted tangential velocity profiles show almost linear 
increases with radial distance out to 10 mm. The exper- 
imental data shows near zero, or even negative tangential 
velocity out to a radius of 5 mm and then a fairly linear 
profile to the tangential velocity peak. This predictive 
behavior is also shown in previous plots. 

The comparisons at an axial distance of 150 mm are 
shown in figure 5.102. Baseline simulations overpredict the 
peak by almost 500 degrees. Pdf simulations slightly 
overpredict the temperature peak. All simulations show a 
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strong drop in temperature towards the combustor axis, while 
temperature measurements show a nearly constant profile. 
Turbulent kinetic energy is best predicted in the pdf simu- 
lations. The baseline simulations overpredict centerline 
turbulent kinetic energy and axial velocity. The jet seeded 
axial velocity measurements are larger than the annulus 
seeded velocity measurements at the same radii. The pdf 
simulations predict axial velocity profiles lying between 
the sets of conditional measurements. Fine grid predictions 
of tangential velocity are very similar to each other, as 
are the two coarse grid simulations. The fine grid simula- 
tions best predict tangential velocity, while the coarse 
grid simulations overpredict tangential velocity by 50 
percent . 

The largest disagreement in the sets of conditional 
velocity measurements occurs for radial velocity, which is 
shown in figure 5.103. The experimental data does not 
benchmark or validate the predictions. Instead, the data 
largely bounds the various velocity predictions. Measure- 
ments of fuel seeded flow show peaks at least twice as high 
as all predictions. Radial velocity measurements taken at 
the 75 and 150 mm axial locations by seeding the fuel jet 
are much higher than measurements found by seeding the 
annulus and external air flows. To best compare tradition- 
ally calculated velocities, velocity measurements should be 
unconditional or averaged. All inlet flows should be 
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simultaneously seeded. 



Chapter 6 

CONCLUDING REMARKS 

A joint pdf solver for species and temperature was 
incorporated and further developed. Velocity, pressure, 
turbulent kinetic energy and dissipation are solved for 
using the SIMPLE algorithm of Patankar. The pdf model uses 
Monte Carlo solution technique to solve for the evolution of 
the pdf. The effects of transport between fluid cells, 
mixing and combustion within fluid cells are treated sequen- 
tially. The composition pdf handles the turbulence combus- 
tion interaction due to varying species and temperature. 

This eliminates the need to correct reaction rate terms due 
to turbulence effects, as is typically done in current 
combustion models. Turbulence causes large local variation 
in species and temperature. As reaction rates are non- 
linear functions of species and temperature, reaction rates 
determined from mean species and temperature are physically 
incorrect. Pdf modeling allows the proper incorporation of 
reaction rates without corrections. 

The pdf method used in this work employs a constant 
number of particles in each fluid cell. An Eulerian frame- 
work is to solve for the evolution of the pdf. The solution 
of particle transfers, mixing and combustion gives a finite 
number representation of the true pdf. This corresponds 
nicely to measurement techniques which also use a finite 
number of data points to ascertain pdfs. 
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Various improvements were made to the pdf model. 

First, variable C p and proper treatment of thermodynamic 
properties was incorporated into the model. The original 
model assumed a single specific heat for all species and 
temperatures . Temperature rather than enthalpy was used to 
keep track of energy. The improved thermodynamic model uses 
data in the form of a two temperature range fifth order 
polynomial to calculate specific heat and enthalpy. This 
way of handling thermodynamic data has been widely used in 
equilibrium and finite rate single cell or one dimensional 
combustion calculations. The method is accurate over a wide 
temperature range and is very highly regarded in the thermo- 
dynamic field. An improved Newton-Raphson species solver 
was incorporated into the composite pdf for solving finite 
chemical reaction rates . 

The validation of the improved hybrid pdf model in- 
volves the comparison of pdf simulations, baseline calcula- 
tions, and experimental data. The baseline combustion model 
uses the SIMPLE algorithm to solve for velocity, pressure, 
species and enthalpy. The baseline model also includes the 
effect of variable specific heat. Originally it was desired 
to calculate spray combustion, so the baseline combustion 
model was evaluated for these flows. The liquid phase fuel 
spray was solved in an Lagrangian framework. The interac- 
tion between the liquid phase and gaseous phase was handled 
by including drag effects on liquid drops and source terms 


in the gas phase for the evaporating spray. The baseline 
spray predictions for an evaporating methanol case were very 
good. Acceptable results were also obtained for the iso- 
thermal swirling air in a spray combustor by including some 
turbulence modeling changes. The baseline model also gave 
qualitatively good results for the same combustor with 
heptane spray combustion. Turbulence model changes did not 
need to be done for the combusting case, as turbulence 
levels generated by combustion were higher. The region 
containing spray was not particularly well calculated de- 
spite various changes. The largest source of error in the 
spray calculation was caused by an overly quick dissipation 
of the gaseous velocity components surrounding the spray. 

The spray occurs in a high shear region where the gaseous 
flow is severely skewed to the computational grid. Velocity 
predictions further from the spray atomizer were much bet- 
ter. As this deficiency in calculating spray combustion 
•would not be improved by turbulence combustion models, 
further testing was confined to single phase flows. It was 
decided not to pursue pdf modeling of combusting sprays. 

There are a few misconceptions about pdf modeling in 
the CFD community. Pdf modeling does not predict emissions 
if they are not included in the reaction mechanisms. Also, 
pdf modeling will not make up for def icientcies in turbu- 
lence modeling. The hybrid pdf model which was used here 
used the k-£ turbulence model to predict mixing. The k-£ 
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model is known to be deficient for high swirl and turbulence 
transition. 

The hybrid pdf was first tested for an enclosed diffu- 
sion combustor. A set of co-annular jets supplied methane 
and air to a dump combustor. Simulations predicted an outer 
recirculation zone. The pdf simulation predicted a recircu- 
lation zone length almost 25% longer. As the recirculation 
zone was hot, and the central fuel core was cool, this 
provided an continuing combustion source. Both simulations 
predicted substantial overlap of fuel and oxygen concentra- 
tion profiles in agreement with measurements. Hybrid pdf 
predictions were most often superior to results obtained 
with the baseline method, especially along the combustor 
centerline. The hybrid pdf produced superior predictions of 
oxygen, carbon dioxide and water concentration. Fuel was 
significantly overpredicted by both hybrid and baseline 
simulations at three axial stations. Both simulations did 
an excellent job of predicting the first temperature pro- 
file. Temperatures on and near the combustor axis were much 
better predicted by the hybrid pdf simulation for the second 
and final temperature profile. The methane combustion 
simulations were beginning to show larger temperature dif- 
ferences with axial distance. More temperature measurements 
would have helped to better evaluate the combustion models. 
Also, measured velocities were not available to validate the 
calculations. Temperature plots of the pdf particles at 
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selected cells showed substantial numbers of partially 
reacted particles co-existing with very hot, fully reacted 
particles. 

The second case involved a swirling combustor flow 
characteristic of an actual gas turbine combustor. The high 
swirl component causes a central recirculation zone. Inlet 
fuel and air flows were at fairly high velocity at a severe 
angle to the computational grid for this combustor. This is 
where the baseline method was identified as being deficient 
for spray calculations. The fuel jet was very coarsely 
modeled. Combustor predictions of velocity and species were 
more diffuse than measurements. The hybrid pdf results for 
the second case did not display as much radial displacement 
of fuel as was measured. Thus, the pdf modeling did not 
make up for the deficiency identified in the baseline com- 
bustor spray calculations. It is likely that a pdf calcula- 
tion done in a Lagrangian framework would prove to be less 
diffusive, and improve flow predictions. Substantial grid 
modification to better model the fuel jet would also help 
modeling this flow. Velocity and turbulence predictions for 
this type of flow would probably be improved with the imple- 
mentation of higher-order numerical scheme. However, such 
fixes are not known for their robustness. The hybrid pdf 
simulation did predict a slower rate of combustion more in 
agreement with experimental temperature and species data. 

The third case involved premixed combustion in a 
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simple dump type combustor. This tested whether pdf model- 
ing could be used in cases were the fuel and oxidant are 
intimately mixed. The hybrid pdf results were generally 
superior to the baseline simulation. The baseline calcula- 
tion best predicted velocity and temperature profiles for 
the first two axial stations. The baseline calculation 
predicted lower temperatures than measured for substantial 
portions of the flow at succeeding axial stations. The 
baseline simulation of combustion seemed overly limited by 
the eddy-breakup reaction rate for latter portions of the 
flowfield. Efforts to improve the baseline prediction by 
increasing the reaction rates resulted in significant change 
in the flowfield and overall solution degradation. The pdf 
simulation predicted higher temperature gradients than the 
baseline simulation. This particularly improved predictions 
for the latter portion of the combustor flowfield. A test 
was done to evaluate the molecular mixing submodel in the 
pdf model. A second simpler relax-to mean mixing submodel 
was incorporated. This change didn't significantly affect 
the pdf predictions. The simpler model is more economical 
to use. 

Very low levels of turbulence were measured in the 
experimental apparatus. A slight deficientcy in the turbu- 
lence model is suggested. While the pdf model produced 
improved temperature profiles, it should not be implied that 
the pdf model make up for def icientcies in turbulence model. 
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The fourth case involved a swirling hydrogen-air jet 
flame. The swirl mass flow component was small enough that 
recirculation was not involved. The hybrid pdf solver pro 
duced superior predictions of temperature, axial velocity, 
and turbulent kinetic energy. The pdf simulation and base- 
line coarse grid simulation both predicted higher tangential 
velocity than measured later in the flowfield. A fine grid 
pdf calculation better predicted tangential velocity decay. 
The baseline fine grid calculation better predicted tangen- 
tial velocity decay. 

The pdf simulations correctly predict turbulent phenome- 
na, specifically, large variation in' local species concen- 
tration and temperature. This represents an improvement in 
turbulent-combustion modeling where turbulent combustion has 
been either ignored or various ad-hoc corrections have been 
applied to chemical reaction rates. Current multi-dimen- 
sional CFD modeling generally produces predict too large 
temperatures and superfluous product species concentrations. 
The inclusion of multiple species concentrations and temper- 
atures ala pdf modeling helps correct this deficiency. The 
pdf method usually improves combustor flowfield predictions. 
These cases prove that the pdf method can be used for both 
diffusion flames and premixed flames provided a proper 
reaction model is used. 

All simulations here used a single reaction step. 
Including multiple reaction rates is known to improve pre- 
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dictions of species and temperature for zero or one dimen- 
sional flows, and allows the calculation of additional 
species, specifically, pollutants. The inclusion of multi- 
ple finite-rate reaction steps in current models, while 
possible, involves increasing numerical inaccuracy. Multi- 
ple reaction steps can involve the introduction of highly 
reactive intermediate species. The solution of reaction 
rates at higher temperatures in current models is only accu- 
rate for very small time steps. Errors generated by using 
too high of reaction rates can produce physically impossible 
predictions of temperature and product species concentra- 
tions. Use of finite reaction rate models in current multi- 
dimensional models involves overwriting negative species 
concentration. Such corrections depend a lot on the experi- 
ence of the CFD modeler. The pdf method can easily incorpo- 
rate chemical kinetic rate solvers which have been developed 
to overcome these problems . These solvers can handle the 
largest chemical kinetic reaction schemes that have been 
proposed. The drawback to this is a rather large amount of 
computer resources and time to perform calculations. Pdf 
modeling may be used to improve the solution accuracy of 
combustor predictions done with simpler combustion models. 


266 


References 


forest 

1. Patankar, S.V., " Numerical Heat Transfer and Fluid Flow", 

Hemisphere Publishing Corp., New York, 1980. 

2. Chen, K-H, Duncan, B., Fricker, D. , Lee, J., and Quealy, A., 
"ALLSPD-3D, Version 1.0", NASA Lewis Research Center, Nov., 1995. 

3. Joshi, D.S., andVanka, S.P., "A Multigrid Calculation Procedure 
for Internal Flows in Complex Geometries", AIAA 90-0442, 1990. 

4. Tafti, D.K., andVanka, S.P., "Hot Gas Environment Around STOVL 
Aircraft in Ground Proximity: Part 2: Numerical Study", AIAA 90- 
2270, 1990. 

5. Vanka, S.P., and Kraz inski, J.L., "Efficient Computational Tool 
for Ramjet Combustor Research", J. Propulsion , Vol. 5, pp. 431-437, 
1989. 

6. Radhakrishnan, K. , and Bittker, D., "LSENS, A General Chemical 
Kinetics and Sensitivity Analysis Code for Homogeneous Gas-Phase 
Reactions", NASA RP-1328, 1329, and 1330, 1994. 

7. Jones, W.P, and Launder, B.E., "The Prediction of Laminarization 

with a two-equation model of turbulence", Int. J. Heat — Mass 

Transfer , Vol 15, pp. 301-314, 1972. 

8. Correa, S.M., and Shyy, W. , "Computational Models and Methods 
for Continuous Gaseous Turbulent Combustion" , Prog. Energy Combust. 
Sci. , Vol. 13, pp. 249-292, 1987. 

9. Leonard, B.P., "A Stable and Accurate Convective Modeling 
Procedure Based on Quadratic Upstream Interpolation" , Compute — Me th . 
Aool. Mech. Engng. , Vol. 19, pp. 59-98, 1979. 

10. Faeth, G.M. , "Current Status of Droplet and Liquid Combustion", 
Proa. Energy Combust. Sci. , Vol. 3, pp. 191-224, 1977. 

11. Law, C.K., "Recent advances in Droplet Vaporization and 
Combustion", Prog. Energy Combust. Sci. , Vol. 8, pp. 171-201, 1982. 

12. Faeth, G.M. , "Vaporization and Combustion of Sprays", Prog^ 
Energy Combust. Sci. , Vol. 9, pp. 1-76, 1983. 

13. Sirignano, W.A. , "Fuel Droplet Vaporization and Spray 
Combustion Theory" , Prog . Energy Combu st . Sci . , Vol. 9, pp. 291 
322, 1983. 

14. Onuma, Y. , and Ogasawara, M., "Studies on the Structure of a 
Spray Combustion Flame", Fifteenth Symposium (International) — on 
Combustion, The Combustion Institute, Pittsburgh, pp. 453-465, 


267 


1975. 

15. Spaulding, D.B. , "The Combustion of Liquid Fuels", Fourth 
Symposium (International) on Combustion, The Combustion Institute, 
Pittsburgh, pp. 847-864, 1953. 


16 Faeth, G.M., "Mixing, Transport and Combustion in Sprays", 
linn Energy Combust. Sci. . Vol . 13, pp. 293-345, 1987. 

17. Tsai, J.S., and Sterling, A.M. , "The Combustion of Linear 
Droplet Arrays " , Twenty-Third — Symposium ( International ) on 
Combustion , The Combustion Institute, Pittsburgh, pp. 1405 141 , 
1990. 


18. Odgers, J. , Kretschmer, D., Norster, E.R., Lefebvre, 
Fletcher, R.S., " The Design and De velopment — of — Gas 
Combustors " , Northern Research and Engineering Corp., 
Massachusetts, 1980. 


A.H. , and 
Turbine 
Woburn, 


19. Tong, A.Y., and Sirignano., W.A., 
Droplet Vaporization with internal 
Formulation and Approximate Solution , 
Vol. 10, pp. 253-278, 1986. 


"Multicomponent Transient 
Recirculation: Integral 

Numerical Heat Transfer, 


20 Raju, M.S., and Sirignano, W.A. , "Multicomponent Spray 
Computations in a Modified Centerbody Combustor", J. of Propulsion , 
Vol. 6, pp. 97-105, 1990. 


21. Reid, R.C., Prausnitz, J.M., and Sherwood, T.K., " Properties of 
Liouids and Gases ", McGraw-Hill, 1977. 

22. Kumar, R. K. , "Ignition of Hydrogen_Oxygen_Diluent Mixtures 
Adjacent to a Hot, Nonreactive Surface," Combustion and Flame, Vol. 
75, pp. 197-215, 1989. 

23. Kim, S.-W., "Numerical Investigation of Chemical Reaction- 
turbulence Interaction in Compressible Shear Layers," Combustio n 
and Flame , Vol. 101, pp. 197-208, 1995. 

24. Dryer, F.L., and Glassman, I., "High Temperature Oxidation of 
Co and CH 4 ", Fourteenth Symposium (Inter national) — on Combustion , 
The Combustion Institute, Pittsburgh, pp. 987-1003, 1972. 


"Simplified 


25. Westbrook, C. K. , and Dryer, F. L. 

for the Oxidation of Hydrocarbon Fuels m 

and Technology, Vol. 27, pp. 31-43, 1981 


Mechanisms 

Combustion 


Reaction 
Flames" , 


Science 


26 Coffee T P., "On Simplified Reaction Mechanisms by Oxidation 
of Hydrocarbon Fuels in Flames by C. K. Westbrook and F. T. Dryer", 
rombust. Sci. and Tech., Vol. 43, pp. 333-339, 1985. 



268 


27. Coffee, T. P. , Kotlar, A. J., Miller, M. S., "The Overall 
Reaction Concept in Premixed, Laminar, Steady-State Flames. I. 
stoi chiometries Combustion and Flame , Vol. 54, pp. 155-169, 1983. 

28. Nikjooy, M. , So, R.M.C., "On the Modelling of Non-Reactive and 
Reactive Combustor Flows", NASA CR 4041, 1987. 

29.Srivatsa, S. K. , "Computations of Soot and Nox Emissions from 
Gas Turbine Combustors", NASA CR-167930, 1982. 

30. Hautman, D.J., Dryer, F.L., Schug, K.P., and Glassman, I., A 
Multi-step Overall Kinetic Mechanism for the Oxidation of 
Hydrocarbons", Combust. Sci. and Tech. , Vol. 25, pp. 219-235, 1981. 

31. Abdalla, A.Y. , Bradley, D. , Chin, S.B. , and Lam, C., ''Global 
Reaction Schemes for Hydrocarbon Oxidation" , Oxidation 
Communications , Vol. 4, pp. 113-130, 1983. 

32. Duterque, J., Borghi, R., and Ticht insky, H. , "Study of Quasi- 
Global Schemes for Hydrocarbon Combustion" , Comb. S ci. and Tech^, 
Vol. 26, pp. 1-15, 1981. 

33. Westbrook , C.K., and Dryer, F.L., "Chemical Kinetic Modeling 
of Hydrocarbon Combustion", Prog. Energy Comb ust. Sci., Vol. 10, 
pp. 1-57, 1984 

34. Paczko, G. , Lefdal, P.M., and Peters, N. , "Reduced Reaction 
Schemes for Methane, Methanol, and Propane Flames ", Twenty- first 
Symposium (International) on Combustion , The Combustion Institute, 
Pittsburgh, pp. 739-748, 1986. 

35. Seshadri, K. , Peters, N. , "Asymptotic Structure and Extinction 
of Methane-Air Diffusion Flames", Combustion and Flame , Vol. 73, 
pp. 23-44, 1988. 

36. Kiehne, T.M., Matthews, R.D., and Wilson, D.E., "An Eight-Step 
Kinetics Mechanism for High Temperature Propane Flames", Combust. 
Sci . and Tech. , Vol. 54, pp. 1-23, 1987. 

37. Jones, W.P., and Lindstedt, R.P., "Global Reaction Schemes for 
Hydrocarbon Combustion", Combustion and Flame , Vol. 73, pp. 233- 
249, 1988. 

38. Bilger, R.W., and Starner, S.H., "On Reduced Mechanisms for 
Methane-Air Combustion in Nonpremixed Flames", Combustion — and 
Flame , Vol. 80, pp. 135-149, 1990 . 

39. D. B. Spaulding, "Mixing and Chemical Reaction in Steady, 
Confined Turbulent Flames" , Thirteenth Symposium (International) on 
Combustion , The Combustion Institute, Pittsburgh, pp. 649-657, 
1970. 


Ill 1 


269 


40. D. B. Spaulding, "Concentration Fluctuations in a Round 
Turbulent Free Jet," Chem. Ena. Sci , vol 26, pp. 95-107, 1971. 

41. Mason, H. B. , and Spaulding, D. B. , "Prediction of reaction 
rates in turbulent pre-mixed boundary layer flows" , Combustion 
Institute European Symposium , pp. 601-606, 1973. 

42. Spaulding, D. B. , "Mathematical Models of Turbulent Flames; A 
Review," Comb. Sci. and Tech. , Vol 13, pp. 3-25, 1976. 

43. Magnussen, B. F., and Hjertager, B. H., "On Mathematical 
Modeling of Turbulent Combustion With Special Emphasis^ On Soot 
Formation and Combustion, " Sixteenth Symposi um (Inte rnational) on 
Combustion , The Combustion Institute, Pittsburgh, pp. 719-729, 
1976. 

44 Srinivasan, R., et. al., "Aerothermal Modeling Program Phase I 
Final Report," NASA CR-168243, GARRETT 21-4742-1, 1983. 

45. Libby, P. A., and Williams, F. A., " Turbulent Reacting Flows, " 
Springer-Verlag, 1980. 

46. Richardson, J.M. , Howard, H.C., and Smith, R.W. , "The Relation 
Between Sampling Tube Measurements and Concentration Fluctuations 

in a Turbulent Gas Jet", Fourth Symposium (International) — on 

Combustion , The Combustion Institute, Pittsburgh, pp. 814-817, 
1953. 

47. Jones, W. P. , "Models for Turbulent Flows with Variable 
Density and Combustion, " in Predic t ion Methods for Turbulent Flows 
(Edited by Kollman, W. ) , pp. 379-422, Hemisphere Publishing Corp., 
1980. 

48. Tsuji, H. and Yamoka, L. "Structure Analysis of Counterflow 
Diffusion Flames in the Forward Stagnation Region of a Porous 
Cylinder, " Thirteenth Symposium (International) on Combustion , The 
Combustion Institute, Pittsburgh, pp. 723-731, 1971. 

49. Bilger, R. W. , "Reaction Rates in Diffusion Flames", Combustion 
and Flame , Vol. 30, pp. 277-284, 1977. 

50. Jeng, S.-M., "An Investigation of Structure and Radiation 
Properties of Turbulent Buoyant Diffusion Flames, Ph.D. Thesis, 
Pennsylvania State University, 1984. 

51. Spaulding, D. B., " GENMIX: A General Computer Program for Two- 
dimensional Parabolic Phenomena ", Pergamon Press, 1977. 

52. Liew, S. K., Bray, K. N. C., and Moss, J. B.,_ "A flamelet model 
of Turbulent Non-premixed Combustion", Comb. Sci. — and Tech .. , Vol. 
27, pp. 69-73, 1981. 


53. Liew, S . k. , and Bray, K.N.C., "A Stretched Laminar Flamelet 
Model of Turbulent Nonpremixed Combustion", Combustion and Flame , 
Vol. 56, pp. 199-213, 1984. 

54. Drake, M.C., and Blint, R.J., "Structure of Laminar Opposed- 
flow Diffusion Flames With CO/H 2 /N 2 Fuel", Combust. Sci. and Tech. , 
Vol. 61, pp. 187-224, 1988. 

55. Borghi, R. , "Turbulent Combustion Modelling", Rrog^. — Energy 
Combust . Sci. , Vol. 14, pp. 245-292, 1988. 

56. Peters, N. , "Laminar diffusion flamelet models in nonpremixed 
turbulent combustion". Prog. Ene rgy Combust. Sci., Vol. 10, pp. 
319-339, 1984. 


57. Janika, J. , and Ko liman, W. , "A two-variable formalism for the 
treatment of chemical reactions in turbulent H 2 -air diffusion 
flames. Seventeenth Symposium (International,} — on Combustion , Tne 
Combustion Institute, Pittsburgh, pp. 421-431, 1979. 

58. Dixon-Lewis, G. , Goldsworthy, F. A., and Greenbur gh , J . B . , 
"Flame structure and flame reaction kinetics, IX. Calculation o 
properties of multi-radical premixed flames", Proc. R. Soc . Lon d^, 
Vol. 346, 261-278, 1975. 

59 Janika, J., and Kollman, W. , "The Calculation of Mean Radical 
Concentrations in Turbulent Diffusion Flames", Combustion and 
Flame , Vol. 44, pp. 319-336, 1982. 

60. Kent, J.H., and Bilger, R.W. , "Turbulent Diffusion Flames", 
Fourteenth Symposium (Internation al) on Combustion, The Combustion 
Institute, Pittsburgh, pp. 615-625, 1974. 


61. Correa, 
"Prediction 
diffusion 
Combustion, 


S.M., Drake, M.C., Pitz, R.W. , and Shyy, W., 
and Measurement of a non-equilibrium . turbulent 

flame". Twentieth Symposium (In ternational) on 

The Combustion Institute, Pittsburgh, pp. 337-343, 


1984. 


62. Correa, S.M., "A model for non-premixed turbulent combustion of 
CO/H 2 jets", Archvm. Combust. , Vol. 5, 223-242, 1985. 


63. Pope, S.B., and Correa, S.M., 
equilibrium Turbulent Diffusion 
(International) on Combus tion, 
Pittsburgh, pp. 1341-1348, 1986. 

64. Mongia, H.C., "Combustion 
Applications and Future Direction 


Joint PDF Calculations of a Non- 
Flame", Twenty- fir st Symposium 
The Combustion Institute, 


Modeling in Design Process: 
", AIAA 94-0466, 1994. 


65. Kuo, K.K., 
York, 1986. 


Principles of Combustion ", John Wiley & Sons, New 


271 


66. Williams, F.A., " Combustion Theory ", Benjamin/Cummings 
Publishing Company, Inc., Menlo Park, 1985. 

67. Bartok, W. , and Sarofim, A.F., "Fossil — Fuel — Combustion , — A 
Source Book", John Wiley & Sons, New York, 1991. 

68. Chigier, N., " Energy. Combustion and Environment", McGraw-Hill, 
New York, 1981. 


69. Oran, E.S., and Boris, J.P, " Numerical Si mulatio n of Reactiv e 
Flow" , Elsevier, New York, 1987. 


70. Anand, M.S., Pope, S.B., and Mongia, H.C., 
Turbulent Recirculating Flows", in Turbulent 
Lecture Notes in Engineering, Springer— Ver lag, 
693, 1989. 


"A PDF Method for 
Reactive Flows , 
Vol. 40, pp. 672- 


71 Curl R L , "Dispersed Phase Mixing, I, Theory and Effects of 
Simple Reactors " , A.I.Ch.E.J. , vol. 9, p. 175-181, 1963. 


72. Janika, J. , Kolbe, 
Transport Equation for 
Turbulent Scalar Fields, " 
47-66, 1979. 


W., and Ko liman, W. “Closure of the 
the Probability Density Function of 
J. Non-Ecruilib- Thermodyn. , Vol. 4, pp. 


73. Dopazo, C., 
jet: Centerline 


Pdf approach for a turbulent axisymmetric heated 
evolution". Physics of Fluids , Vol. 18, 397-404, 


1975. 


74. Pope, S. B., "Pdf Methods for Turbulent Reactive Flows", Prog^ 
Energy Combust. Sci. , Vol. 11, pp. 119-192, 1985. 


75. Hsu, A. T., and Chen, J. Y. , "A Continuous Mixing Model for Pdf 
Simulations and its Application to Combusting Shear Flows , iq 
symposium on Turbulent Shear Flows , Technical University of Munich, 
p. 22-4, Sept. 9-11, 1991. 


76 Chen, J.Y., and Kollman, W. , "Segregation Parameters and Pair 
Exchange Mixing Models . for Turbulent Nonpremixed Flames", Twent y^ 

Third Symposium (International) on Combustion, The Combustion 

Institute, Pittsburgh, pp. 751-757, 1990. 


77. Pope, S. B. , "An Improved Turbulent Mixing Model", Comb j — Sci. 
Tech. . Vol. 28, pp. 131-135, 1982. 


78. Norris, A. T. , and Hsu, A., 
Closure Methods in the Modeling of 
TM 106614, 1994. 


"Comparison of PDF and Moment 
Turbulent Reacting Flows", NASA 


79. Masri, A. R. , Dibble, R. W., and Barlow, R. S., "The Structure 
of turbulent, pilot stabilized Flames of CY{JCO/n 2 /^ fuel 
mixtures", Combustion and Flame , Vol. 88, pp. 239-264, 1992. 


80. Villermaux, J., "Micromixing phenomena in stirred reactors", 
Encyclopedia of Fluid Mechanics , pp. 707-768, Gulf Pub. Co., West 
Orange , N J , 1986. 

81. Yamazaki, and Ichigawa, A., Inti. Chem. Eng. , Vol. 10, pp. 471- 
478, 1970. 

82. Correa, S.M., and Braaten, M.E., "Parallel Simulations of 
Partially Stirred Methane Combustion", Combustion and Flame , Vol. 
94, pp. 469-486, 1993. 

83. Correa, S.M., "A Direct Comparison of Pair-Exchange and I EM 
Models in Premixed Combustion", Combustion and Flame , Vol. 103, pp. 
194-206, 1995. 

84. J.Y. Chen, and Kollman, W. , "Mixing Models for Turbulent Flows 
with Exothermic Reactions", Turbulent Shear Flows 7 , Springer- 
Verlag, Berlin Heidelberg, pp. 277-292, 1991. 

85. Chen, J. Y., Kollman, W. , and Dibble, R. W. , "Pdf modeling of 
Turbulent Nonpremixed Methane Jet Flames, " Combust. Sci. and Tech. , 
Vol. 64, pp. 315-346, 1989. 

86. Maas, U.A. , and Pope, S. B., "Simplifying chemical kinetics: 
Intrinsic Low-Dimensional Manifolds in Composition Space", 
Combustion and Flame , Vol. 88, pp 239-264, 1992. 

87. Chen, J.Y. , and Kollman, W. , "Chemical Models for Pdf Modeling 
of Hydrogen-Air Nonpremixed Turbulent Flames", Combustion — and 
Flame . Vol. 79, pp 75-99, 1990. 

88. Hsieh, K.-C., Shuen, J.-S., Tsai, Y.-L. Peter, and Yu, S.-T., 
"RPLUS2D/3D USER'S MANUAL," Sverdrup Technology, Inc., NASA Lewis 
Research Center, Oct. 22, 1990. 

89. Yoon, S., and Jameson, A. "An LU-SSOR Scheme for the Euler and 
Navier-Stokes Equations," AIAA Paper 87-0600, 1987. 

90. Burros, M. C., and Kurkov, A. P., "Analytical and Experimental 
Study of Supersonic Combustion of Hydrogen in a Vitiated Air 
Stream," TMX-2828, NASA, 1973. 

91. Rogers, R. C., and Chinitz, W. , "Using a Global Hydrogen-Air 
Combustion Model in Turbulent Reacting Flow Calculations", AIAA J . , 
Vol. 21, pp. 586-592, 1983. 

92. Strahle, W.C., and Lekoudis, S.G., "Evaluation of Data on 
Simple Turbulent Reacting Flows", AFSOR-TR-85-0880 , 1985. 

93. Kline, S.J., Morkovin, M.J., and Moffat, "Computation of 
Turbulent Boundary Layers", 1968 AFSQR-IFP-Stanfor d Conference, 
1969. 


273 


and Whitelaw, J.H. 


94. Baker, R.J., Hutchinson, P., Khalil, E.E. . 

"Measurements of Three Velocity Components in a Model Furnace With 

and Without Combustion", Fifteenth Symposium — ( International ) Qft . 

Combustion , The Combustion Institute, Pittsburgh, pp. 553 559, 
1974. 


95. Lockwood, F.C., El-Mahallawy , F.M., and Spaulding, D.B., "An 
Experimental and Theoretical Investigation of Turbulent Mixing in 
a Cylindrical Furnace", Combustion an d Flame, Vol. 23, pp. 283 y , 
1974. 


96. McDannel, M.D. , Peterson, P.R. , and Samuelson, G.S., _ Species 
concentration and Temperature Measurements in a Lean, Premixed Flow 
Stabilized by a Reverse Jet", Combust. Sci. and Tech , Vol. 28, pp. 
211-224, 1982. 


97. El Banhawy, Y., Sivasegaram, S., and Whitelaw, J.H. , "Premixed, 
Turbulent Combustion of a. Sudden- Expans ion Flow" , Combustion a nd 
Flame, Vol. 50, pp. 153-165, 1977. 


98. Smith, G.D., Giel, T.V., and Catalano, C.G.^, 
Reactive Recirculating Jet Mixing in a Combustor", 
pp. 270-276, 1983. 


"Measurements of 
AIAA J. , Vol. 21, 


99. Lightman, A. J. , Richmond, 
P.D. , Roquemore, W.M. , Bradley 
C.M., "Velocity Measurements 
AIAA-80-1544, 1980. 


R.D., Krishnamurthy, L., Magill, 
, R.P., Stutrud, J.S., and Reeves, 
in a Bluff-Body Diffusion Flame", 


100. Magill, P.D., Lightman, A.J., Orr, 
Roquemore, W.M. , "Simultaneous Velocity 
in a Bluff-Body Combustor", AIAA-82-0833 


C.E., Bradley, R.P., and 
and Emission Measurements 
, 1982. 


101. Schefer, R.W., Namazian, M., and Kelly, J. , 
Measurements in a Turbulent Nonpremixed Bluff-Body 
Flame", AIAA-87-1349 , 1987. 


"Velocity 

Stabilized 


102. Sivasegaram, S. and Whitelaw, J.H., "Oscillations in 
Axisymmetric Dump Combustors", Combust. Sci. and Tech ^., Vol. 52, 
pp . 413-426, 1987. 

103. Sivasegaram, S., Thompson, B.E., and Whitelaw, J.H., "Acoustic 
Characterization Relevant to Gas Turbine Augmentors , J^ 
Propulsion, Vol. 5, pp. 109-115, 1989. 


104. Charles, R., "Detailed Data Set: Velocity and Temperature 
Measurements in the Axisymmetric Can Combustor for a Parametric 
Variation in Inlet Conditions", University of California, Irvine, 
UCI-ARTR-87-6 , 1988. 


274 


105. Pan, J.C., "Laser-Diagnostic Studies of Confined Turbulent 
Premixed Flames Stabilized by Conical Bluff Bodies", Ph.D. Thesis, 
University of Dayton, 1991. 

106. Heitor, M.V., and White law, J.H., "Velocity, Temperature, and 
Species Characteristics of the Flow in a Gas-Turbine Combustor", 
Combustion and Flame , Vol. 64, pp. 1-32, 1986. 

107. Bicen, A.F., Tse, D.G.N, and Whitelaw, J.H., "Combustion 
Characteristics of a Model Can-Type Combustor" , Combustion and 
Flame , Vol. 80, pp. 111-125, 1990. 

108. Bulzan, D.L. "Structure of a Swirl-Stabilized Combusting 
Spray", Journal of Propulsion and Power , Vol. 11, pp. 1093-1102, 
1995. 

109. Srivasan, R. , "Aero thermal modeling program Phase I final 
report", NASA CR- 1682 43, 1983. 

110. Sloan, D.G., Smith, P.J., and Smoot, L.D. , "Modeling of swirl 
in turbulent flow systems". Proa. Energy Combust. Sci. , Vol. 12, 
163-250, 1986. 

111. Shih, T-H, Zhu, J. , and Lumley, J.L., "Modeling of wall- 
bounded complex flows and free shear flows," NASA TM 106513, 1994. 

112. Me Done 11 , V.G., Adachi, M. , and Samuel sen, G.S., "Structure of 
Reacting and Non-Reacting Swirling Air-Assisted Sprays", Combust. 
Sci. and Tech. , Vol. 82, pp. 225-248, 1992. 

113. El Banhawy, Y., and Whitelaw, J.H., "Calculation of the Flow 
Properties of a Confined Kerosene- Spray Flame", AIAA. J . , Vol. 18, 
pp. 1503-1510, 1980. 

114. Edwards, C.F., and Rudoff, R.C., "Structure of a Swirl- 
Stabilized Spray Flame by Imaging, Laser Doppler Velocimetry, and 
Phase Doppler Anemometry" , Twenty-Third Symposium (International) 
on Combustion , The Combustion Institute, Pittsburgh, pp. 1353-1359 , 
1990. 

115. Bachelo, W.D., Brena de la Rosa, A., and Rudoff, R.C., 
Diagnostics Development for Spray Characterization in Complex 
Turbulent Flows", AIAA 88-0236, 1988. 

116. Zhu, J.Y., Rudoff, R.C., Bachalo, E.J. and Bachalo, W.D., 
"Number Density and Mass Flux Measurements Using the Phase Doppler 
Particle Analyzer in Reacting and Non-Reacting Flows", AIAA 93- 
0361, 1993. 

117. Presser, C., Gupta, A.K., and Semerjian, H.G. , "Aerodynamic 
Characteristics of Swirling Spray Flames: Pressure-Jet Atomizer", 
Combust, and Flame, Vol. 92, pp. 25-44, 1993. 



275 


118. Ghaffarpour, M. , and Chehroudi, B., "Experiments on Spray 

Combustion in a Gas Turbine Model Combustor", Combust . Sci . und 

Tech. , Vol . , 92, pp. 173-200, 1993. 

119. Lewis, M. H., and Smoot, L. D. "Turbulent Gaseous Combustion 
Part I: Local Species Concentration Measurements", Combustion and 
Flame , Vol. 42, pp. 183-196, 1981. 

120. Smith, P. J., and Smoot, L. D. "Turbulent Gaseous Combustion 
Part II: Theory and Evaluation for Local Properties", Combus tion 
and Flame , Vol. 42, pp. 277-285, 1981. 

121. Jones, W. P., and Wilhelmi, J., "Velocity, Temperature and 
Composition Measurements in a Confined Swirl Driven Recirculating 
Flow", Combust. Sci. and Tech. , Vol. 63, pp. 13-31, 1989. 

122. Jones, W. P., Tober, A., "Velocity, Composition and Tempera- 
ture Fields in an Axisymmetric Model Combustor", Intern ational 
Symposium on Applications of Laser Anemometry to Fluid Mechanics, 
Lisbon, Portugal, p. 3.9, 1988. 

123. Gould, R. D., Stevenson, W. H. , and Thompson, H. D. , "Turbu- 
lence Characteristics of an Axisymmetric Reacting Flow" , NASA CR 
4110, Feb., 1988. 

124. Magnussen, B.F., and Hjertager, B.H., "On Mathematical 
Modeling of Turbulent Combustion with Special Emphasis on Soot 

Formation and Combustion" , Sixteenth Symposium (Int ernational) on 

Combustion , The Combustion Institute, Pittsburgh, pp. 719-729, 
1978. 

125. Magnussen, B.F., "The Eddy Dissipation Concept for _ Turbulent 
Combustion Modelling. Its Physical and Practical Implications", 
Division of Thermodynamics, University of Trondheim, Norway, 1990. 

126. Gran, I. R., Melaaen, M. C., and Magnussen, F., "Numerical 
Simulation of Fluid Flow and Combustion in Gas Turbine Combustors", 
ASME 93-GT-166, International Gas Turbine and Aeroengine Congress 
and Exposition, 1993. 

127. Takahashi, F., Vangsness, M.D., Belovich, V.M. , "Conditional 
LDV Measurements in Swirling and Non-Swirling Coaxial Turbulent Air 
Jets for Model Validation," AIAA 92-0580, 1992. 

128. Durbin, M.D., Vangsness, M.D., Ballal, D.R. , and Katta, "Study 
of Flame Stability in a Step Swirl Combustor," ASME 95-GT-lll, 
1995. 

129. Takahashi, F., Vangsness, M.D., Durbin, M.D., and Schmoll, 
W.J., "Structure of Turbulent Hydrogen Jet Diffusion Flames," ASME 
HTD-Vol. 317-2, Proceedings of the ASME Heat Transfer Division , pp. 
183-193, 1995. 


Appendix 1 

GAS PHASE SOURCE TERMS 


276 


Equation 


Axial momentum 


Radial momentum 


P-9 


dx 
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Oxygen mass fraction 
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Turbulent kinetic energy 
Dissipation 
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Turbulent parameters 
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